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This  report  documents  computational  algorithms  -and- 
developed  computer  programs  for  modeling  two-dimensional  (2- 
D)  and  three-dimensional  (3-D)  seepage  under  dams  and 
groundwater  flow  in  aquifers.  A  method  for  automatically 
generating  flow  nets  for  2-D  applications  and  Cray  YMP, 
scientific  visualization,  and  numerics  for  3-D  problems  are 
emphasized.  Specifically,  techniques  typically  used  by 
aerospace  engineers  are  applied  to  flow  through  porous 
media.  This  research  and  development  was  done  and  this 
report  was  written  at  the  US  Army  Engineer  Waterways 
Experiment  Station  (WES) ,  Information  Technology  Laboratory 
(ITL) ,  Interdisciplinary  Research  Group,  Computer-Aided 
Engineering  Division  (CAED) ,  by  Dr.  Fred  T.  Tracy  in  partial 
fulfillment  of  the  requirements  for  the  degree  of  Doctor  of 
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of  the  exit  point. 

Porosity  or  ratio  of  the  volume  of  voids  to 
the  total  volume. 


N 

{N} 

[N] 

Ne 

Ne 

Nf 


Normal  coordinate  of  a  tangent-normal 
coordinate  system. 

Shape  function  vector. 

Shape  function  matrix. 

Number  of  equipotential  drops. 

Number  of  elements. 

Number  of  flow  paths. 

Shape  function  for  the  ith  node  of  an 
element . 
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Symbol 


Description 


Ni 


Shape  function  for  the  first  node  of  an 
element . 


N2 

Shape  function 
element. 

{N), 

Shape  function 
a  face  of  a 

[N], 

Shape  function 
a  face  of  a 

N, 

Shape  function 
element . 

for  the  second  node  of  an 

vector  for  the  four  nodes  of 
brick  element. 

matrix  for  the  four  nodes  of 
brick  element. 

for  the  ninth  node  of  an 


[P] 

0 

% 

'I'e 

't'E 

'I'j 

'I't. 

’t'TJ 


Matrix  of  partial  derivatives  of  shape 
functions . 

Total  head  or  potential. 

Second  Euler  angle. 

Complex  potential. 

Stream  function. 

Value  of  stream  function  at  point  D. 
Third  Euler  angle. 

Value  of  stream  function  at  point  E. 

Value  of  stream  function  at  point  I. 

Value  of  stream  function  at  point  J. 

Value  of  partial  of  stream  function  with 
respect  to  T  at  pointy  I. 

Value  of  partial  of  stream  function  with 
respect  to  T  at  point  J. 

Total  amount  of  stream  function. 


Constant  value  of  stream  function. 


Symbol 


Description 


’l'2 

q 

Q 

Q' 

{Q} 

Qe 

(QIg 

{Q)u 

Qi 

02 

(Q), 

<0)8 

{Q>9 

r 

{r} 

<^)a 


Constant  value  of  stream  function. 

Flux  density  for  the  partially  penetrating 
well  problem. 

Quantity  of  flow. 

Five-component  variable  for  Euler  Equations 
in  strong  conservative  form. 

Nodal  flow  vector. 

Flow  specified  at  an  external  node. 

Global  flow  vector. 

Global  flow  vector  computed  before  the 
first  iteration  for  unconfined  flow 
problems . 

Nodal  flow  vector  computed  from  the 
unmodified  stiffness  matrix  and  the 
nodal  heads. 

Flow  for  node  1. 

Flow  for  node  2. 

Element  flow  vector  for  the  four  nodes  of  a 
face  of  a  brick  element. 

Element  flow  vector,  excluding  the  internal 
node. 

Element  flow  vector,  including  the  internal 
node . 

Radial  distance  from  the  center  of  the 
well . 

Coordinate  vector. 

Coordinate  vector  at  point  A. 

Coordinate  vector  at  point  B. 
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Symbol 

[^]e 

P 

Q 


s 

S, 


t 

T 

tw 

0e 

01 

02 


U 

U 


Description 

Coordinate  matrix  for  an  element. 

Coordinate  vector  at  node  i. 

Radius  of  the  well. 

Value  of  coordinate  vector  where  the 
pressure  head  is  zero. 

Coordinate  matrix  for  the  four  nodes  of  a 
face  of  a  brick  element. 

Density  of  the  soil-water  complex. 

Concentration  of  water  in  the  soil-water 
complex . 

Parameter  that  varies  between  zero  and  one. 

Specific  storage. 

Value  of  s  where  the  pressure  head  is  zero. 

Final  value  of  s  in  the  exit  line 
computation. 

Time. 

Tangent  coordinate  of  a  tangent-normal 
coordinate  system. 

Thickness  of  a  partially  penetrating  well. 

First  Euler  angle. 

Angle  between  an  equipotential  line  and  the 
normal  to  the  boundary  for  material  type 
No.  1. 

Angle  between  an  equipotential  line  and  the 
normal  to  the  boundary  for  material  type 
No.  2  . 

X  component  of  discharge  velocity. 

Energy. 
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Symbol 

{u) 


V 
0 

{V) 

V. 

J 

V 

0 

{v}o 


u . 

I 

w 

X 


Description 

Nodal  displacements  vector. 

Discharge  velocity  at  the  left  flux 
boundary  of  a  finite  volume  cell. 

Discharge  velocity  at  the  right  flux 
boundary  of  a  finite  volume  cell. 

y  component  of  discharge  velocity. 

Velocity  vector. 

Discharge  velocity  vector. 

Discharge  velocity  vector  at  point  A. 

Discharge  velocity  vector  at  point  B. 

Discharge  velocity  vector  at  point  C. 

jth  component  of  discharge  velocity. 

Discharge  velocity  vector  at  node  n. 

Normal  component  of  discharge  velocity. 

Specified  discharge  velocity. 

Discharge  velocity  vector  at  the  average  of 
the  respective  (x,  y,  z)  coordinates  of 
the  corner  nodes  of  an  element. 

jth  component  of  the  actual  velocity  of  the 
water  particles. 

z  component  of  discharge  velocity. 

X  coordinate. 


Xg  X  coordinate  of  the  first  node  for  the 

computation  of  the  exit  point. 

X|^  X  coordinate  of  the  second  node  for  the 

computation  of  the  exit  point. 


XX 


Symbol 

Description 

X  coordinate  of  the  third  node  for  the 
computation  of  the  exit  point. 

{X}e 

X  coordinate  vector  for  an  element. 

X  coordinate  at  node  i  of  an  element. 

X,  y,  or  z,  depending  on  whether  j  =  1,  2, 
or  3 . 

X  coordinate  of  a  point  on  a  boundary. 

(^0'  Vo) 

Location  of  a  well  in  an  aquifer. 

X  coordinate  of  the  first  node  of  an 
element. 

^2 

X  coordinate  of  the  second  node  of  an 
element. 

X  coordinate  vector  for  the  four  nodes  of  a 
face  of  a  brick  element. 

X, 

X  coordinate  of  the  ninth  node  of  an 
element . 

X' 

X  coordinate  in  a  new  coordinate  system. 

($,  ri,  C) 

Coordinates  for  computational  type  space 
where  finite  elements  are  mapped. 

(ii,  C,) 

Coordinates  for  computational  type  space 
where  finite  elements  are  mapped  at  node 

i . 

y 

y  coordinate. 

Ya 

y  coordinate  of  the  first  node  for  the 
computation  of  the  exit  point. 

Yb 

y  coordinate  of  the  second  node  for  the 
computation  of  the  exit  point. 

Yc 

y  coordinate  of  the  third  node  for  the 
computation  of  the  exit  point. 

xxi 

Symbol 


Description 

Value  of  y  coordinate  at  point  D. 

y  coordinate  vector  for  an  element. 

Value  of  y  coordinate  at  point  E. 

y  coordinate  at  node  i  of  an  element. 

y  coordinate  of  a  point  on  a  boundary. 

y  coo^-dinate  vector  for  the  four  nodes  of  a 
face  of  a  brick  element. 

y  coordinate  of  the  ninth  node  of  an 
element . 

y  coordinate  in  a  new  coordinate  system. 

z  coordinate. 

z  coordinate  vector  for  an  element. 

z  coordinate  at  node  i  of  an  element. 

Maximum  z  value  of  the  nodes  of  the  FE 
grid . 

Minimum  z  value  of  the  nodes  of  the  FE 
grid . 

z  value  at  the  top  of  an  aquifer. 

z  coordinate  vector  for  the  four  nodes  of  a 
face  of  a  brick  element. 

z  coordinate  of  the  ninth  node  of  an 
element . 

Normalized  z  for  the  solution  to  a 
partially  penetrating  well. 

Function  used  in  the  solution  to  a 
partially  penetrating  well. 
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CHAPTER  I 


INTRODUCTION 

Introduction 

The  modeling  of  seepage  under  dams  and  groundwater  flow 
in  aquifers  is  of  significant  interest.  This  becomes  even 
more  important  in  our  modern  times  wich  increasing  interest 
in  the  flow  of  pollutants.  The  unsolved  environmental 
issues  regarding  our  hazardous  and  toxic  waste  problems  must 
be  resolved,  and  significant  resources  must  be  placed  on 
this  effort.  Some  military  bases  are  contaminated  with 
hazardous  waste  that  has  entered  the  groundwater  domain.  A 
groundwater  model  that  takes  into  account  contaminant  flow 
is  therefore  critical.  The  state  of  the  art  has  advanced  in 
various  ways  over  the  years  to  achieve  better  and  better 
solutions.  However,  of  unusual  occurrence  is  the  applica¬ 
tion  of  the  tools  that  engineers  in  one  discipline  have 
developed  to  problems  of  other  disciplines.  What  is  said 
is,  "We  don't  do  it  that  way."  Because  of  the  author's 
diverse  background,  a  unique  feature  of  the  work  in  this 
dissertation  is  that  the  tools  developed  by  structural  and 
aerospace  engineers  are  applied  to  a  problem  typically 
addressed  by  others. 

Scope  of  Dissertation 

This  dissertation  concentrates  effort  on  the 
complicated,  real-world  problem  of  seepage  and  groundwater 
flow  in  three  ways: 

1.  The  development  and  application  of  new  and 

innovative  computational  techniques  for  a  more 
effective  solution  procedure. 
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2.  The  application  of  techniques  and  software 
developed  by  structural  and  aerospace  engineers  to 
a  civil  engineering  problem  (technology  transfer) . 

3.  The  development  and  application  of  a  three- 
dimensional  (3-D)  seepage  package  and  groundwater 
model  using  the  latest  grid  generation  and  scien¬ 
tific  visualization  techniques. 

More  detail  will  now  be  given  to  the  different  parts. 


Two-Dimensional  (2-D)  Improvements 

A  2-D  finite  element  method  (FEM)  seepage  package  has 
been  developed  by  the  author  which  has  three  parts  as 
follows : 

1.  Interactive  graphics  grid  generation. 

2.  Steady-state  seepage  analysis  for  both  confined  and 
unconfined  problems,  homogeneous  or  inhomogeneous 
media,  . 'otropic  or  anisotropic  soil,  plane  and 
axisymmetric  flow. 

3.  Interactive  graphics  postprocessor. 

This  package  has  been  distributed  worldwide  and  has  advanced 
the  state  of  the  art.  Further  advancements  in  automatic 
flow  net  generation  for  a  2-D  problem  from  the  perspective 
of  aerospace  engineering  will  first  be  presented.  In  fact, 
the  technique  of  using  the  Cauchy-Riemann  equations  to  gen¬ 
erate  a  structured  orthogonal  grid  is  extended  to  generate  a 
seepage  flow  net.  Several  examples  and  comparison  with 
other  work  are  also  presented. 


Three-Dimensional  Model 

Next,  a  3-D  seepage  package  and  groundwater  model  is 
presented  where  the  primary  hardware  configuration  is  the 
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combination  of  a  Cray  YMP  and  a  Silicon  Graphics  Iris  work¬ 
station.  The  individual  parts  will  now  be  described. 

Grid  Generation 

The  program  EAGLE  (written  predominantly  by  Dr.  Joe 
Thompson,  Mississippi  State  University)  has  extensive  capa¬ 
bility  in  generating  structured  grids  for  finite  volume  flow 
applications.  Many  users  of  FEM  programs  prefer  structured 
grids  (although  not  required)  because  of  their  esthetics  and 
good  numerical  properties  (triangular  and  tetrahedral  ele¬ 
ments  show  bias  at  times) .  Therefore,  EAGLE  was  used  as  the 
grid  generation  tool. 

Conversion  and  Data  Completion 

A  module  to  convert  the  output  from  the  grid  generation 
program  to  the  FEM  seepage  format  is  next  described.  This 
includes  combining  data  from  different  blocks  into  a  single 
FEM  grid,  and  boundary  condition  and  soil  property  data  must 
also  be  supplied.  Finally,  a  bandwidth  minimization  algo¬ 
rithm  must  be  applied  to  the  grid. 

FEM  Seepage  Analysis 

A  3-D  version  of  the  2-D  seepage  analysis  program  is 
next  presented.  The  numerical  problems  in  going  from  2-D  to 
3-D  are  also  described.  Like  the  new  2-D  seepage  program, 
the  capabilities  are: 

No  initial  guess  of  the  free  surface  or  any 
restrictions  on  the  grid  shape  or  numbering  in  the 
vicinity  of  the  free  surface  will  be  required. 
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2.  Layers  with  different  soil  properties  are  allowed. 

3.  The  program  will  terminate  upon  convergence  without 
the  user  having  to  specify  a  specific  number  of 
iterations  or  time  steps. 

Scientific  Visualization 

The  new  3-D  seepage  program  outputs  data  compatible 
with  FAST  (available  from  NASA  Ames  and  operational  on  the 
Iris  Workstation) .  The  techniques  to  properly  format  the 
unconfined  flow  data  are  described.  Since  EAGLE'S  grids  can 
be  plotted  with  FAST,  graphics  capability  for  both  prepro¬ 
cessing  and  scientific  visualization  are  supplied  to  the 
user. 

Test  Problems 

The  developed  3-D  seepage  and  groundwater  model  was 
tested  with  both  theoretically  verifiable  and  real-world 
problems.  These  results  are  presented. 

Previous  Work 

Early  attempts  at  modeling  groundwater  (Meyer  and 
Kleinecke  1968)  assumed  horizontal  flow  in  cells  of  one 
layer  (2.5-D)  with  a  finite  difference  scheme.  However, 
when  the  flow  became  unconfined,  the  problem  became  nonlin¬ 
ear,  and  it  was  difficult  to  ensure  convergence  in  a  steady- 
state  problem.  In  fact,  one  significant  aspect  of  seepage 
is  handling  the  problem  of  unconfined  flow  with  a  free 
surface  through  materials  of  significantly  different  charac¬ 
teristics  (permeabilities) .  One  example  of  this  is  an  earth 
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dam  with  a  relatively  impervious  clay  core  with  a  highly 
pervious  drain  installed  around  it.  The  rest  of  the  dam  is 
composed  of  moderately  porous  material.  Figure  1  shows  an 


example 

with  the 

following  soil 

properties . 

Permeabilities,  ft/min 

Soil 

Material 

^^2 

Angle,  deg 

1 

Rock 

9. 3  (10‘2) 

1.7(10'^) 

140 

2 

Sand 

9.8(10'^) 

2.0(10'^) 

0 

3 

Drain 

9.8(10'^) 

2.0(10‘2) 

0 

4 

Shell 

9.8  (10'^) 

9.8(10'') 

0 

5 

Random 

9.8(10'^) 

9.8(10'^) 

0 

6 

Core 

9.2(10'^) 

2.0(10'^) 

0 

7 

Random 

9.8(10'^) 

9.8(10'^) 

0 

8 

Grout 

9.8(10'^) 

2.0(10'^) 

140 

A  finite  element  analysis  (FEA)  (Wilson  1969, 
Zienkiewicz  1971)  for  2-D  steady-state  seepage  flow  for 
confined  or  unconfined  flow  (Taylor  and  Brown  1967,  Finn 
1967)  was  developed.  However,  the  user  was  required  to 
orient  the  elements  where  the  free  surface  might  occur  in  a 
special  way  and  give  angles  along  which  the  free  surface 
should  move.  This  was  an  early  attempt  at  adaptive  mesh 
methods  but  did  not  work  well  at  times.  Also,  this  code 
could  not  handle  problems  that  became  partially  confined  and 
partially  unconfined  unless  the  user  could  pick  a  priori 
where  the  break  would  be.  Some  improvement  to  the  conver¬ 
gence  problem  was  made  (Neuman  and  Witherspoon  1970) .  A 
different  approach  to  a  finite  difference  solution  for 
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;  -y  /  v  /  /  /  /  /-y  /  /  /  /  /^  /“7  7 '7  /////// 

Figure  1.  Earth  Dam  with  Different  Soil  Properties 

a  transient  problem  of  an  initially  dry  bank  (Desai  1970) 
was  taken,  but  this  solution  was  good  for  only  one  boundary 
condition  and  only  a  rectangular  mesh.  Later,  a  method  of 
solving  the  transient  seepage  problem  using  the  FEM  by 
treating  the  problem  as  a  series  of  steady-state  problems 
(France  1971)  was  developed.  Further  refinements  (Issacs 
and  Mills  1972)  to  the  work  of  France  also  were  made  by 
modifying  elements  crossed  only  by  the  free  surface  rather 
than  "adapting"  an  entire  set. 

At  this  time  the  author  developed  both  2-  and  3-D  FEM 
seepage  programs  (Tracy  1973a,  1973b) .  These  programs  were 
used  in  a  study  of  Lock  &  Dam  26,  Alton,  IL,  with  reasonable 
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results  (Hall,  Tracy,  and  Radhakrishnan  1975).  However,  a 
subsequent  project  failed  where  the  grid  was  produced 
manually. 

It  was  apparent  that  grid  generation  and  interactive 
graphics  tools  were  essential  for  the  successful  application 
of  numerical  techniques  to  large  real-world  problems.  Work 
was  then  begun  on  a  2-D  FEM  grid  generator  and  postprocessor 
(Tracy  1977a,  1977b,  1977c) .  The  grid  generator  had  two  new 
capabilities  as  follows: 

1.  After  the  user  defined  points  and  line  segments, 
the  program  automatically  put  them  together  into 
subregions.  If  they  were  four-sided  subregions,  a 
structured  algebraic  grid  was  generated. 

2.  If  the  subregion  had  an  arbitrary  number  of  sides, 
a  triangular  mesh  was  generated  by  a  simple,  yet 
innovative,  technique  developed  by  the  author. 

The  postprocessor  could  do  the  following  for  2-D  FEM  output: 

1.  Numbers. 

2.  Contours. 

3.  Vectors. 

4.  Displaced  grid  for  structures  problems. 

5.  Isometric. 

6.  Perspective. 

As  PC's  and  engineering  workstations  became  more  powerful, 
Apollo  and  IBM  PC  versions  were  created  (Tracy  1988) .  Work 
was  done  to  plot  a  3-D  FEM  grid  (Tracy  and  Wade  1980) ,  but 
the  work  to  generate,  a  3-D  grid  was  left  to  others.  That 
is,  in  fact,  one  of  the  reasons  why  the  unique  capability  of 
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EAGLE  (Thompson  1987,  Thompson  and  Gatlin  1988a,  1988b, 
1988c)  is  used  in  the  work  documented  by  this  dissertation. 

The  seepage  programs  developed  at  this  point  have  some 
serious  limitations: 

1.  The  scheme  to  do  a  steady-state  problem  by  allowing 
the  transient  problem  to  converge  with  a  constant 
time  increment  is  very  slow. 

2.  When,  as  shown  in  Figure  2,  the  free  surface 
crosses  the  grid,  small  and  therefore  sometimes 
skewed  elements  just  below  the  free  surface  are 
created . 

3.  With  materials  with  significantly  different 
permeabilities,  a  time  step  good  for  one  layer  is 
totally  wrong  for  the  others. 

4.  An  initial  guess  of  the  free  surface  is  required. 

5.  The  elements  where  the  free  surface  would 
potentially  go  have  to  be  quadrilaterals  numbered 
in  a  certain  way,  destroying  one  of  the  major 
features  of  the  FEM. 


Figure  2.  FEM  Grid  with  Free  Surface 


A  new  approach  where  the  elements  above  the  free  surface  are 
given  a  very  low  permeability  (Bathe  and  Khashgoftaar  1979) 
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was  developed  which  alleviated  many  of  these  problems. 

There  is  one  boundary  condition  (at  the  exit  point)  which 
was  approximated  in  the  work  of  Bathe  and  is  very  difficult 
to  handle.  A  modified  version  of  Bathe's  work  with  a 
correct  version  of  exit  point  boundary  conditions  (Tracy 
1983)  was  done  next.  This  program  was  later  put  into  the 
above-mentioned  2-D  seepage  package  (Biedenharn  and  Tracy 
1987)  with  the  pre-  and  postprocessor  programs.  The  last 
improvement  is  work  on  the  addition  of  a  flow  net  option 
(Tracy  and  Radhakrishnan  1989)  which  is  described  in  this 
dissertation  with  a  new  aerospace  engineering  perspective. 
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CHAPTER  II 


GOVERNING  EQUATIONS 

Flow  Ecmations 

The  equation  for  unsteady  unconfined  flow  of  a  com¬ 
pressible  fluid  in  an  initially  dry  compressible  porous 
medium  (Dewiest  1966)  such  as  in  the  earth  dam  shown  in 
Figure  3  with  headwater  level  Hj,{t)  and  tail  water  level 


where 


p  =  mass  density  of  the  fluid 

u,  V,  w  =  discharge  velocities  in  the  x,  y,  and  z 
directions,  respectively 

n  =  porosity  (ratio  of  the  volume  of  voids  to  the 
total  volume) 

t  =  time 

The  actual  velocities  of  the  Ti’jid  particles  are  related  to 
the  discharge  velocities  by 


where 

Uj  =  jth  component  of  the  actual  velocity  of  the  water 
particles 

Vj  =  jth  component  of  discharge  velocity  (u,  v,  or  w) 
Further,  the  concentration  of  the  water  in  the  soil-water 
complex  can  be  measured  by  a  density  g  defined  by 

g  =  np  (2.3) 

Equations  2.2  and  2.3  can  now  be  substituted  into  Equation 
2.1  to  yield  the  conservation  of  mass  equation  (Anderson, 
Tannehill,  and  Fletcher  1984)  used  by  the  aerospace  engi¬ 
neers  as  follows: 

T  "  4^  =  0  (2.4) 

dXj  dt 

where 

X.  =  X,  y,  or  z 

This  equation  can  be  written  in  vector  form  by 
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V  •  (eu)  +  If  =  0 

where  u  is  the  velocity  vector.  In  like  manner. 
Equation  2 . 1  can  be  written 


V  -  (pv)  +  =  0 

where  v  is  the  discharge  velocity  vector.  If  the  soil-water 
complex  is  assumed  incompressible,  p  is  constant.  Also,  a 
more  general  equation  can  be  stated  for  either  confined  or 
unconfined  flow.  The  above  two  steps  yield 


V  •  V  + 


0 


where 

Sg  =  specific  storage 
h  =  total  head  or  potential 

Darcy's  Law  (Maasland  1957)  for  laminar  flow  for  discharge 
velocity  is  given  by 


V  =  -  JC  •  Vh 


(2.5) 


where  K  is  a  3  X  3  permeability  tensor.  So  now  the  equation 
for  unsteady  flow  becomes 


V  ■  (K  •  Vh) 


(2.6) 


Steady-State  Solution 

Because  this  work  solves  the  difficult  problem  of  3-D 
unconfined  flow  with  no  restriction  on  the  grid,  automatic 
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flow  net  generation,  and  building  a  working  and  tested  tool 
for  use  by  industry  and  federal  installations,  it  concen¬ 
trates  on  a  steady-state  solution.  If  a  steady-state  solu¬ 
tion  is  approached  from  solving  a  nonlinear  equation,  then 
the  following  equation  is  used- 

V  ■  (K  •  Vh)  =  0  (2.7) 

Finally,  in  the  special  case  where  a  homogeneous,  isotropic 
medium  is  modeled,  the  total  head  satisfies  Laplace's 
Equation , 

=  0  (2.8) 

Equation  2.8  will  be  used  to  generate  a  flow  net  for  2-D 
problems . 

! 
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CHAPTER  III 


COMPUTER  GENERATED  FLOW  NETS 

Introduction 

The  graphical  construction  of  flow  nets  by  hand  to 
compute  the  quantity  of  flow,  exit  gradient,  etc.  is  a  stan¬ 
dard  engineering  tool  of  soils  engineers.  However,  these 
are  extremely  tedious  to  construct  by  hand  because  equipo- 
tential  lines  and  flow  lines  must  be  drawn  in  such  a  way 
that  curvilinear  squares  result.  One  significant  aspect  of 
this  research  effort  is  that  numerical  grid  generation  tech¬ 
niques  of  aerospace  engineers  used  to  generate  an  orthogonal 
grid  (Thompson,  Warsi,  and  Mastin  1985)  can  be  extended  to 
construct  a  flow  net  for  various  boundary  conditions  using 
the  Cauchy-Riemann  Equations  (Crowder  and  McCuskey  1964). 
This  chapter  shows  how  the  FEM  has  been  successfully  applied 
to  generate  flow  nets  with  emphasis  also  given  to  differ¬ 
ences  in  approach  from  previous  work.  The  major  advantage 
of  the  techniques  described  in  this  chapter  is  that  they 
improve  the  quality  of  the  resulting  flow  nets. 


Governing  Equations  and  Basic  Approach 
The  total  head  or  potential  0  for  a  homogeneous, 
isotropic  medium  for  2-D,  steady-state  flow  (as  already 
shown  by  Equation  2.8)  satisfies  Laplace's  equation  as 
follows : 


^  ^  =  0 
dx^  dy^ 


(3.1) 

The  stream  function  i|i  also  satisfies  Laplace's  Equation, 


^  ^  =  0 
dx^  dy^ 

Therefore,  a  complex  potential  i  exists  as  follows; 


<()  =  4>  +  iijf 


(3.2) 


Since  0  and  41  are  conjugate  harmonic  functions,  the  Cauchy- 
Riemann  equations  now  hold. 


dx 


dy 


or  = 


'I'y 


04)  _  _  0l|f 

dy  dx 


or  4>y  = 


(3.3) 


It  should  be  noted  here  that  the  stream  function  is  often 
defined  as  a  velocity-type  term.  However,  in  this  work  it 
is  a  gradient-type  term.  That  is,  \|;  is  defined  by 


as  compared  to 


4)  =  4)^dy  -  (\>ydx 


udy 


vdx 
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A  property  of  such  functions  is  that  constant  lines  of 
0  are  orthogonal  to  constant  lines  of  ij;.  The  flow  net  con¬ 
sists  of  0  =  constant  lines  and  iji  =  constant  lines 
constructed  in  such  a  way  that  the  resulting  picture  con¬ 
sists  of  curvilinear  squares.  The  concept  of  automatically 
generating  the  flow  net  is  fairly  straightforward  and 
involves  the  following  three  steps: 

1.  Perform  a  normal  FEM  solution  determining  the  total 
head  h  (same  as  potential  0)  at  each  node  and  the 
quantity  of  flow,  Q,  passing  through  the  system. 
Also  compute  the  shape  factor,  f,  from 

Q  =  k{h^  -  f  (3.4) 

where  h^  is  the  upstream  head,  h^^  is  the  downstream 
head,  k  is  the  permeability,  Q  is  the  flow,  and  f 
is  the  shape  factor. 

2 .  Determine  the  boundary  conditions  for  the  stream 
function  and  perform  a  second  FEM  solution  to 
obtain  values  of  i|f  at  each  node. 

3.  Contour  the  two  sets  of  data  to  construct  the  flow 
net.  The  intervals  for  each  are  determined  using 
the  shape  factor  which,  by  definition,  is 


where  is  the  number  of  equipotential  drops,  and 
Nf  is  the  number  of  flow  paths. 

Earlier  work  (Christian  1980a,  1980b,  1983;  Aalto  1984; 

Christian  1987)  determined  the  boundary  values  for  the 

stream  function  solution  (step  2)  numerically,  whereas  in 

this  work  a  more  fundamental  technique  is  used.  Here,  the 
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Cauchy-Riemann  equations  are  used  to  determine  the  correct 
boundary  conditions. 


Modified  Cauchv-Riemann  Conditions 
A  modified  form  of  Equation  3.3  will  be  applied  to 
determine  the  boundary  conditions  for  the  stream-function 
computation  (step  2) .  On  the  boundary,  a  tanqent-normal 
(T-N)  coordinate  system  is  used  as  illustrated  in  Figure  4 
It  is  created  by  first  constructing  a  parallel  coordinate 
system  x'  -  y'  at  the 

y 


I - - 

Figure  4.  Tangent-Normal  System 

boundary  point  (x^,  y^)  and  then  rotating  the  prime  system 
an  angle  A.  The  local  and  global  coordinate  systems  are 
related  by: 


X 


X 


y'  =  y  -  Vo 

and  the  tangent-normal  coordinate  system  is  related  to  the 
x'  -  y'  system  by: 


T]  f  cos  A  sin  A\(x 


N)  \-sin  A  cos  A  [v^ 


where  A  is  the  angle  between  the  two  axes.  Therefore, 


>t\  cos 


A  sin  aV4>, 


I4)J  1,-sin  A  cos  AM 


(3.4) 


'IfrI  (  cos  A  sin  aV’J' 


lij;  J  v-sin  A  cos 


(3.5) 


where 


dT  dT 


dN  dN 


Applying  Equation  3.3  to  Equation  3.5  yields 


'I'rl  (  cos  A  sin  A)f~*t>y 


\-sin  A  cos  A/[  <))^ 
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Collecting  terms  gives 


1 

■e 

_ 

cos  A  sin  A 
-sin  A  cos  A 


<t> 

\^y/ 


(3.6) 


Since  the  right-hand  side  of  Equations  3.4  and  3.6  are  the 
same,  the  left-hand  sides  can  be  equated  to  obtain 


4>t  = 

<l>Af  =  -’l»r 


(3.7) 


Equation  3.7  can  now  be  used  to  solve  for  the  boundary  con¬ 
ditions  for  the  stream  function. 

Examples 

Several  examples  will  now  be  presented  to  illustrate 
the  procedure. 


Four-Sided  Regions 

Two  problems  showing  examples  of  the  basic  four-sided 
region  are  now  given.  They  are  flow  under  a  weir  and  two 
versions  of  flow  in  a  partially  penetrating  slot. 

Confined  Flow  under  a  Weir 

This  simple  problem  (Figure  5)  clearly  demonstrates  the 
concepts  presented  in  this  chapter.  Note  that  there  are 
four  boundaries,  with  constant  head  specified  on  two  lines 
(segments  AFE  and  BCD)  and  a  no-flow  (impervious)  condition 
specified  on  the  other  two  lines  (segments  AB  and  ED) . 
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Figure  5.  Weir  Problem 

On  the  0  =  and  0  =  boundaries  (AFE  and  BCD) 

4>t  =  0 

since  and  Hj  are  constants.  But  from  Equation  3.7  we  get 

<t»T  = 

So  the  new  boundary  condition  is 

=  0 

On  the  impervious  boundaries  no  flow  enters,  so  the 
normal  component  of  velocity,  v^,  is  zero.  Thus,  using 
Darcy's  Law  for  a  homogeneous,  isotropic  medium, 

=  0 

it  follows  that 

Applying  Equation  3.7  to  the  above  equation  yields 
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<t‘N  =  'Vr  =  - 


or 

=  0  (3-8) 

Now 

d\\f=}]fj.dT+y\i^dN  (3.9) 

Substituting  Equation  3.8  into  the  above  equation  and  noting 
that  dN  =  0  on  the  boundary  gives 

di|r  =  0 
or 

ilf  =  constant  (3.10) 

The  total  amount  of  stream  function  can  be  shown  to  be 

4'cotai  =  -f  =  -  *d) 

Thus  for  the  impervious  boundaries  (using  Equation  3.10), 
set 

vjr  =  =  0  on  AB 

'I'  =  4^2  =  total  on  ED 

(In  the  actual  computer  program,  a  constant  is  added  to 
and  il»2  so  their  values  will  be  higher  than  the  maximum  y 
value  of  the  grid  to  keep  the  FE  program  from  solving  an 
unconfined  flow  problem.) 

Notice  that  the  boundary  conditions  are  exactly 
reversed.  The  impervious  boundary  in  the  first  problem  is 
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replaced  as  a  constant  head  boundary,  and  a  constant  head 
boundary  has  been  replaced  as  a  no-flow  or  impervious 
boundary . 

A  wide  variety  of  problems  involving  four-sided  regions 
can  be  handled  by  the  concept  demonstrated  here.  The  first 
example  is  given  in  Figure  6  where  a  weir  similar  to  the  one 
of  Figure  5,  but  with  sheet  piles  added,  is  shown  (BC  =  AJ  = 
40  ft,  DE  =  GH  =  10  ft,  CD  =  IJ  =  30  ft,  and  FG  =  40  ft) . 
Figure  7  shows  the  computer  generated  flow  net  for  this 
problem.  The  sheet  piles  remain  part  of  the  one  continuous 
impervious  boundary  (DEFGHI  in  Figure  6) .  Points  D  and  F 
have  the  same  (x,  y)  coordinates  as  do  points  G  and  I. 


Figure  6.  Weir  Problem  with  Sheet  Piles 

Partially  Penetrating  Slot 

The  second  example  is  confined  flow  through  a  partially 
penetrating  slot.  Both  the  case  of  specified  head  and 
specified  constant  discharge  velocity  at  the  slot  are 
considered . 
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Head  specified  at  the  slot 

Figure  8  shows  that  head  or  potential  is  specified  at 
two  boundaries  (at  the  headwater  level  on  BC  and  at  the  slot 
along  DEF) ,  and  the  other  two  boundaries  are  impervious 
(along  the  center  line  and  base  FAB  and  the  top  CD) .  Since 
this  problem  is  topologically  identical  to  the  weir  problem, 
no  further  work  is  needed  to  solve  it. 


H - 7-  J - 7 - 7“ 

Figure  8. 


7  / - 7 - 7  7  V 


c 


H: 


~7-  J  r  — 7 - 7"'  J  -V - o 

Partially  Penetrating  Slot 


Constant  discharge  velocity 
specified  at  the  slot 

This  version  of  the  slot  problem  consists  of  constant 

discharge  velocity  specified  on  one  of  the  four  boundaries, 

a  constant  head  on  one  boundary,  and  impervious  conditions 
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on  the  other  two  (Figure  9) .  Line  segments  EAB  and  DC  are 


Figure  9.  Specified  Discharge  Velocity 


impervious,  line  segment  BC  has  constant  head  specified, 
and  line  segment  DE  has  a  specified  discharge  velocity,  V^. 
On  the  0  =  H,  boundary,  use  as  before, 

=  0 

On  the  boundary  DE, 

or 


Equation  3.7  can  again  be  applied  to  obtain 


so 


Starting  again  with  Equation  3.9, 

diir  =  j  dT  *  ^  dN 
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and  observing  that  dN  =  0  on  the  boundary  again  gives 


Integrating, 


'I' 


+  c 
k 


where  c  is  a  constant  to  be  evaluated.  In  this  particular 
problem  y  coincides  with  T  along  DE,  so 


Let 


c 


then 


at  y  =  Ye 


c 


and 


^  ^  iy  -  Ye)  * 

Now,  on  the  impervious  boundaries,  apply  the  = 
constant  boundary  condition  as  follows: 

vlf  =  i|f£  on  EAB 
ijr  =  ij/p  on  CD 

where  \|ip  is  computed  by 
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'I'd  =  ^  (Yd  -  Ye^  ^  I'd 

The  boundary  conditions  for  the  stream-function  solution 
are  now  completely  defined. 

Dupuit's  Problem 

Dupuit's  problem  deals  with  unconfined  flow  in  an  earth 
dam  with  vertical  sides  (Figure  10) .  Line  segment  AB  is 
impervious,  line  segments  AF  and  BC  have  constant  head  spec¬ 
ified,  and  line  segment  CD  has  the  boundary  condition 

4>  =  y 

The  position  of  the  free  surface  (FD)  must  be  determined 
from  the  FEM  solution.  Once  determined,  line  segment  FD 
becomes  a  flow  line  and  is  treated  exactly  like  an  impervi¬ 
ous  boundary.  Also,  the  region  above  the  free  surface 


(triangular  region  FDE)  is  not  used  for  the  second  solution. 
Rather,  a  new  grid  with  the  phreatic  surface  being  a  new 
boundary  is  used.  We  will  now  determine  the  i|i  boundary 
conditions . 

On  the  0  =  and  0  =  H2  boundaries  use,  as  before, 

=  0 

On  the  impervious  boundaries  apply  the  41  =  constant  boundary 
condition  as  follows: 

4;  ==  on  AB 

4f  -  =  4/1  +  4ftotai  on  FD 

On  the  boundary  CD, 

<t>  =  y 

This,  however,  is  insufficient  information  to  determine  the 
new  boundary  conditions.  Therefore,  the  normal  component  of 
discharge  velocity  is  first  computed  for  each  node  on  the 
boundary  CD.  It  is  assumed  that  points  C  and  D  are  node 
points,  and  there  are  typically  intermediate  node  points 
such  as  nodes  I  and  J  in  Figure  10  as  well.  Then  for  each 
node , 

=  -'I'r 
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It  is  further  assumed  that  (and  therefore  i|f^)  varies 
linearly  between  node  points.  So  between  two  node  points  I 
and  J  having  \|i  values  of  41,  and  il»j  and  i|t^  values  of  \|i^|  and 
is  approximated  as 


'^T 


Integrating  as  before  gives 

Solving  for  4ij  gives 

^  trj)  +  ’I'l 

The  first  set  of  node  points  starts  with  point  I  corre¬ 
sponding  to  point  D  (Figure  10) .  Thus, 

'I'j  =  '^D 

Also,  using 

=  ^”  =  0 

(since  is  zero  all  along  the  flow  line  FD) ,  begin  with 

ilf„  -  0 

Point  D,  being  the  exit  point,  causes  significant  numerical 
problems  (to  be  discussed  in  Chapter  IV) .  With  this  start 
and  the  values  of  computed  for  all  the  remaining  nodes  on 


28 


the  surface  of  seepage  from  the  FEM  solution,  the  nodes 
along  DC  can  now  be  processed  consecutively  until  the  tail- 
water  point  C  is  processed.  The  boundary  conditions  are  now 
fully  determined  for  the  stream-function  calculation.  Note 
that  the  only  restriction  on  is  that  it  be  large  enough 
to  prevent  the  FEM  analysis  program  from  solving  an  uncon¬ 
fined  flow  problem  while  solving  for  the  stream  function. 

The  above  formulation  is  not  restricted  to  Dupuit's 
Problem  but  can  also  be  applied  to  a  wide  variety  of  quadri¬ 
lateral-type  earth  dams.  The  only  restriction  is  that  the 
problem  should  have  the  five  basic  boundaries  of  imper”^lous 
base,  specified  headwater  and  tailwater,  free  surface,  and 
surface  of  seepage.  Figure  11  shows  the  results  for 
Dupuit's  Problem  where  by  Figure  10,  AB  =  AF  =  100  ft,  and 
BC  =  20  ft.  Figure  12  shows  the  computer  generated  flow  net 
for  an  earth  dam. 


Figure  11.  Results  for  Dupuit's  Problem 
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Anisotropic  Soil 

Flow  nets  for  an  anisotropic  soil  can  be  produced 
equally  well  for  one  layer  with  this  technique.  The  number 
of  flow  lines  and  equipotential  lines  is  computed  the  same 
way  as  with  the  isotropic  case.  However,  the  resulting  plot 
will  now  have  curvilinear  quadrilaterals  instead  of  curvi¬ 
linear  squares.  To  get  curvilinear  squares,  the  geometry 
and  permeability  must  be  transformed  in  the  traditional  way 
as  follows: 

x"  =  X 


y 


N 


y 
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k 


where  and  are  the  horizontal  and  vertical  perme¬ 
abilities,  respectively.  If  the  principal  permeabilities 
are  not  coincident  with  the  x-y  axis,  x  and  y  must  represent 
a  rotated  coordinate  system  in  the  equations,  except  that 
now  kj^  and  k^  are  more  properly  rendered  k^  and  k2. 

Figure  13  shows  an  anisotropic  problem  consisting  of  a  weir 


Figure  13.  Weir  on  Anisotropic  Soil 


with  a  sheet  pile  and  a  sloping  impervious  bottom.  Note 
that  the  principal  permeabilities  are  not  parallel  to  the 
x-y  axis.  Figure  14  shows  the  flow  net  if  the  soil  is  as¬ 
sumed  isotropic,  and  Figure  15  shows  the  results  for  the 
anisotropic  case. 
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Figure  14.  Isotropic  Results 


Multilayered  Problems 

Multilayered  problems  where  all  the  flow  lines  origi¬ 
nate  or  end  in  the  same  material  type  can  also  be  solved. 
To  understand  how,  first  consider  the  two-layered  system 
shown  in  Figure  16.  The  procedure  is  done  the  same  as  the 
single-layered  case  except  that  in  the  second  solution  for 
flow  lines  the  permeabilities  are  modified.  The  reason  is 
that  when  an  equipotential  or  flow  line  crosses  a  boundary 
between  materials  of  different  permeabilities,  it  is  bent. 
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Figure  16.  Flow  across  Boundary 


However,  since  both  regions  are  isotropic,  the  flow  lines 
must  remain  perpendicular  to  the  eguipotential  lines  on  both 
sides  of  the  boundary.  Flow  lines  are  bent  as  follows 
(Freeze  and  Cherry  1979)  (Figure  16) : 


tan  01  _ 
tan  02  k2 


(3.11) 


In  a  similar  manner,  since  eguipotential  lines  are  perpen¬ 
dicular  to  flow  lines,  eguipotential  lines  are  bent  as 
follows : 


^ 

tan  ( -  Oj) 
or 

cot  Cj  _  k^ 
cot  a2  k^ 

Finally,  taking  the  reciprocal  of  both  sides  gives 
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tan  _  k2 
tan  a 2 


(3.12) 


When  the  second  FEM  solution  is  completed,  the  program 
treats  the  flow  lines  as  equipotential  lines.  Therefore, 
the  flow  lines  are  bent  according  to  Equation  3.12  instead 
of  Equation  3.11.  To  compensate  for  this,  the  ratios  of 
permeability  must  be  reversed.  This  can  be  accomplished, 
for  instance,  by  reversing  k,  and  k^.  For  a  multilayered 
problem,  however,  it  is  best  to  replace  the  k's  with  their 
reciprocals,  respectively.  Thus,  the  ratio  of  k's  between 
material  types  1  and  2  becomes 


This  enables  the  procedure  to  work  for  any  number  of  layers. 

Figure  17  shows  a  three-layered  problem  with  the  divid¬ 
ing  line  between  the  layers  being  equidistant,  and  Figure  18 
shows  the  computer  generated  results.  Since  the  problem  has 
isotropic  soil  layers,  the  equipotentials  and  flow  lines 
remain  perpendicular  to  each  other.  However,  in  one  layer 
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layers  there  will  be  curvilinear  rectangles. 

Axisymroetr ic  Case 

The  flow  net  as  traditionally  defined  (a  plot  of  equal 
potential  drops  and  equal  flow  paths  constructed  in  such  a 
way  as  to  produce  curvilinear  squares)  does  not  exist  for  an 
axisymmetric  problem  (Figure  19  shows  the  result  for  a  fully 
penetrating  well  which  has  a  well  radius  of  1  ft,  a  radius 
of  influence  of  51  ft,  a  depth  of  20  ft,  and  a  permeability 
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of  0.1  ft/min.).  Notice  that  tall,  thin  rectangles  become 
short  and  fat  proceeding  from  the  well  outward.  No  addition 
or  subtraction  of  flow  lines  can  change  this  situation.  The 
technique  outlined  in  this  chapter  can  be  used  to  produce  a 
plot  showing  mutually  perpendicular  equipotential  and  flow 
lines,  but  the  number  of  each  is  arbitrary. 

Comparison  v/ith  Other  Work 
Two  problems  presented  by  Christian  were  dealt  with 
using  the  techniques  presented  in  this  chapter  and  the 
results  compared.  The  first  problem  is  confined  flow  in  an 
anisotropic  medium  as  shown  in  Figure  20.  Figure  20  also 
shows  Christian's  results  with  the  author's  results  in 
Figure  21.  Note  that  the  results  are  the  same.  A  more 
d'fficult  problem  is  unconfined  flow  through  a  two-zoned 
earth  dam  with  kg  =  2k^.  The  problem  and  Christian's 
results  are  shown  in  Figure  22,  and  the  author's  results  are 
shown  in  Figure  23.  With  curvilinear  squares  in  Region  A, 
curvilinear  rectangles  of  size  2  to  1  should  occur  in 
Region  B.  Note  that  the  results  by  the  author  are  much 
closer  to  this  result.  Also,  due  to  numerical  complexities 
at  the  exit  point  that  the  author  avoids,  an  extra  wiggle 
occurs  in  Christian's  work  near  the  exit  point.  Clearly, 
the  more  fundamental  approach  generates  a  superior  product. 
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't>  CONTOUR  INTERVAL  -  1  000 

a  MAXIMUM  .  12.000 
O  MINIMUM  .  0 

V  CONTOUR  INTERVAL  -  2.0000 

O  MAXIMUM  .  15.655 
O.  MINIMUM  .  .0 


Figure  20.  Confined  Flow  in  an  Anisotropic  Medium 


Figure  21.  Author's  Results 
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<p  CONTOUR  INTERVAL  -  10  000 
a  MAXIMUM  -  85.000 
O  MINIMUM  .  0 

y  CONTOUR  INTERVAL  -  10.000 
O  MAXIMUM  .  41  826 
O  MINIMUM  -  .0 
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CHAPTER  IV 

THREE-DIMENSIONAL  NUMERICS 

Introduction 

The  goal  of  this  work  is  to  allow  a  completely  general 
3-D  FE  grid  for  both  confined  and  unconfined  flow  problems. 
Further,  the  user  is  not  to  be  required  to  make  an  initial 
guess  of  free  surface  areas  or  be  concerned  about  the  number 
of  iterations  required  for  convergence.  The  advantage  of 
this  approach  is  that  a  totally  structured  grid,  partially 
structured  and  partially  unstructured  grid,  or  totally 
unstructured  grid  can  be  used.  However,  going  from  2-  to 
3-D  is  significant.  For  example,  in  a  2-D  unconfined  flow 
problem  there  is  typically  one  exit  point  to  cause  computa¬ 
tional  problems.  In  a  3-D  problem  there  is  an  exit  line 
instead  of  an  exit  point,  and  there  can  be  any  number  of 
them.  One  example  is  an  aquifer  or  cofferdam  with  several 
wells.  The  numerics  invo''ved  in  the  3-D  solution  will  now 
be  discussed. 

Finite  Volume  versus  Finite  Element  Comparison 
An  unusual  approach  will  be  used  to  introduce  the  FE 
formulation.  Typically,  a  variational  approach  is  used  to 
derive  the  governing  FEM  equations.  This  works  very  well 
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and  is  at  times  physically  based  such  as  in  structures  prob¬ 
lems  where  the  principle  of  virtual  work  in  variational  form 
can  be  used.  What  will  be  done  here,  however,  is  to  first 
compare  finite  volume  and  FEM  by  considering  one¬ 
dimensional  (1-D),  steady-state  seepage  flow  in  the  region 
shown  in  Figure  24  with  three  finite  volume  cells  having 
dimensions  of  Ax  by  Ay  by  1. 


\/\/\/\/\/  \  X 


# 

# 

• 

i-1 

• 

I 

i  +  1 

/  \  \ 

/  \  /  \ 

/  \  /  \ 

Figure  24. 

One-Dimensional 

Seepage  Flow 

First,  putting 

Equation  2.5  back 

into  Equation  2.7 

gives 


V •  V  =  0  (4.1) 

Integrating  finite  volume  i  over  its  volume  gives 

f  V  ■  V  dV  =  0 

J  V 

Using  the  Divergence  Theorem  produces 
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J  V  ■  ds  =  0 

Since  flow  is  horizontal,  only  the  two  vertical  boundaries 
have  flux.  So  the  equation  can  be  approximated  by 

u  1  Ay  -  u._i  Ay  =  0  (4.2) 

2  2 

The  1-D  version  of  Equation  2.5  for  homogeneous,  isotropic 
flow  is 


dx 

Using  a  linear  version  of  Richardson's  Interpolation  and 
gathering  information  on  either  side  of  the  boundary 
produces  the  standard  approximation  at  the  boundaries  of 


“'4 


^i*l  ~  ^2 

Ax 

Ax 


The  above  series  of  equations  can  be  substituted  in 
Equation  4.2  to  obtain 

-  khll-lJllAy  ^  T-^L-^Ay  =  0 

Ax  Ax 

When  terms  are  collected,  one  obtains 


('1,.,  -  2h,.  0  (4.3) 

The  term  in  parenthesis  is  the  familiar  central  difference 
formulation  for  the  second  derivative  in  computational 


space . 


2 


Now  the  same  problem  can  be  considered  by  using  the 
direct  stiffness  approach  of  the  FEM  (Figure  25) .  The  three 
cells  are  now  three  finite  elements  with  eight  nodes.  Since 
the  flow  is  horizontal,  nodes  1  and  5,  2  and  6,  3  and  7,  and 
4  and  8,  respectively,  have  the  same  head.  Therefore, 
nodes  5  through  8  are  considered  inactive.  A  structures 


problem  will  relate  forces  and  displacements  by  use  of 
Hook's  Law  via  a  stiffness  matrix  as  follows: 

{F}  =  [^r]  {u} 

where 

{F}  ==  nodal  forces 

[K]  =  stiffness  matrix 

{u}  =  are  the  nodal  displacements 

In  an  analogous  manner,  in  seepage  the  nodal  flows  (Q)  are 
linked  to  nodal  heads  as  follows: 
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{Q]  =  m  {h} 


Let  us  now  look  at  element  1.  x  varies  inside  the 
element  as  follows: 


X 


■-^1 

2  ^ 


X, 

2  2 


=  +  N2X2 


where 

-1  ^  ^  1 

The  vector  form  is 

x={N]^{x]^  (4.4) 

where  {N}  is  a  vector  of  shape  functions  and  {x)^  are  the 
nodal  X  coordinates  for  an  element.  Since  nodes  1  and  2  are 
the  only  active  nodes  for  element  1,  let  h  vary  the  same  as 
X  in  Equation  4.4.  Therefore, 

h  =  {N}’^{h}^  (4.5) 

where  {h}^  are  the  nodal  heads  for  an  element.  Thus,  the 
isoparametric  element  formulation  (Cook  1981)  has  been  used. 

The  nodal  flows  are  simply  the  flux  (volume  of  water 
per  unit  time)  crossing  at  that  node,  so  in  this  case, 

Q  =  uAy 

But 

u  =  -ki 
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where  the  gradient  i  is  in  this  case 


From  the  definitions, 


i  = 

dx  ® 


dx  di 


[P]  = 


d{N} 


we  get 


i  =  J-Hp]  ih}, 


Our  goal  here  is  to  obtain  a  matrix  [B]  so  that 


i  =  [B]  {h}, 


(4.6) 


We  see  from  these  equations  that 


[B]  =  J-UP] 


(4.7) 


In  our  1-D  example 


-r  1  /  ,  A;v 

-  -  — 


[P]  =  ^  [-1  1] 
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[B]  =  ^[-1  1] 

Ax 


The  discharge  velocity  can  now  be  computed  as 


u  =  -k[B]  {h}, 


(4.8) 


which  becomes  a  simple  constant  in  the  1-D  problem  as 
follows : 


u  =  (h,  -  h^) 

Ax 

Our  goal  is  to  determine  the  element  stiffness  matrix 
[K]g  relating  the  nodal  Q's  to  the  nodal  h's.  The  sign 
convention  used  for  flow  entering  the  element  is  positive 
and  for  flow  exiting  the  element  is  negative.  Thus, 

Ca  =  -  h,) 


'0,  (  1 

icAy 


So  from  the  definition 


(4.9) 


we  get 
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kAy 

Ax 


This  can  be  written 


[K]  ^  - 


[k]  ^(-1  1) 

Ax 


AxAy 


The  above  equation  can  also  be  written 


[K]^=  [B]^[k]  [B]  AxAy 

So  our  final  goal  is  reached,  which  is  to  have  an  intuitive 
basis  for  the  general  expression  of  the  element  stiffness 
matrix,  which  is 


[iC]^  =  [S]^[ic]  [B]  dV 


(4.10) 


where  now  [k]  is  the  matrix  version  of  the  general  perme¬ 
ability  tensor. 


General  Finite  Element  Method  Formulation 
We  are  now  ready  to  discuss  the  general  FEM  equations 
and  techniques  for  the  seepage  problem. 


Element  Formulation 

The  different  aspects  of  the  individual  element  formu¬ 
lation  will  now  be  discussed. 
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Element  Shape 

To  conform  to  the  output  from  EAGLE,  as  well  as  add  a 
quadratic  variation  in  the  interpolation  function,  a  nine- 
node  element  is  used.  As  shown  in  Figure  26,  this  becomes 
an  eight-node  "brick"  element  with  an  internal  node  having 
coordinate  values 


Figure  26.  Finite  Element  Type 


Unstructured  grids  are  accomplished  by  nodes  collapsing  to 
form  prisms,  tetrahedrons,  etc.  as  shown  in  Figure  27. 


Isoparaunetric  Element 

The  general  versions  of  Equations  4.4  and  4.5  will  now 
be  given.  First,  the  general  finite  element  of  Figure  26  is 
mapped  to  a  2-  by  2-  by  2 -unit  cube  in  ($,  rj ,  C) 
computational  type  space,  where 

-1  ^  ^  ^  1 

•  -1  ^  T1  ^  1 

-1  ^  C  1 

The  values  of  $,  r) ,  and  C  at  the  corner  nodes  are  either  -1 
or  1,  and  (0,  0,  0)  at  node  9.  The  ith  interpolation  func¬ 
tion  Nj  for  values  of  i  from  1  to  8  becomes 
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(1  +  (1  +  CiO  (4.12) 

o 

where  (i-,  r}.,  Cj)  are  the  ($,  rj ,  C)  coordinates  at  point  i. 
The  interpolation  function  at  point  9  is 

Wg  =  (1  -  e)  (1  -  h")  (1  -  C")  (4.13) 

Now  h  is  interpolated  by 

8 

h  =  ^9^9  (4.14) 

i=l 


where 


(4.15) 


SO  that  the  matrix  version  again  becomes  as  in  Equation  4.5, 


Because  {r},,  where 


h  = 


{r} 


y 


is  defined  at  the  average  point  of  the  element,  {r}'  at  node 
9  is  zero.  So, 


{r} 


1  =1 


A  matrix  form  for  this  equation  is 
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{0}^  {0}n[ 

{r}  =  {0}^  {0}^  {y}« 

,{0}^  {0}^ 

In  these  equations, 


Also,  {0)  is  a  zero  vector.  Let  the  above  equation  for  {r} 
be  written  as 

{r}  =  [N]  [r]^ 

where 

'{N}'^  {0}^  {0}^ 

[N]  =  {0}^  {N}^  {0}^  [r]  ^  =  {yle 

,{o}^  {0}’’  [^^K, 

Stiffness  Matrix  Formulation 

The  element  stiffness  matrix  as  given  in  Equation  4,10 
is  now  derived.  First,  the  general  equation  for  the  gradi¬ 
ent  is 
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{i}  =Vh  = 

So  from  Equation  4.6, 

[B]  = 

Now  we  will  use  the  operators 


d\ 

(  ^  \ 

dx 

di 

d 

dy 

= 

_d_ 

dtl 

_d_ 

d 

.  3z, 

\  acj 

to  determine  [B].  First,  the  Jacobian  Matrix  [J]  is 
by 

ij)  =  -If 

d(i,  ti,  C) 

This  can  also  be  written 


[J]  =  V'lr}"' 

=  V{[N]  [r]J^ 

Expanding, 

[J]  {N}^{y}^  {N}^{z}^) 

-  {y}^  {z}^) 


The  matrix. 


[P]  =  V'{N]^ 

has  only  i,  r] ,  and  C  variables,  so  it  can  be  easily 
uated.  Tnarefore, 


(4.16) 


defined 


(4.17) 

eval- 
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[J]  =  [P]  ({x}^  {y}^  {z}J 


(4.18) 


Next,  using  the  chain  rule  of  differentiation, 

V'  =  [J]V 

For  elements  that  are  not  skewed  or  ill-conditioned,  [J]  is 
well  defined.  Therefore, 

V  =  [cJ]  -^v' 

[B]  then  becomes 

[B]  =  [J]^^V{N}^ 

Using  Eguation  4.17  we  get 

[P]  =  [cJj-HP]  (4.19) 

[B]  is  now  completely  defined  in  terms  of  constants  and 
computational  coordinates. 

The  element  stiffness  matrix  as  given  in  Eguation  4.10 
can  now  be  evaluated  by  integrating  in  computational  space 
as  follows: 

[iCle  =  f_[  f_\  [P]  IJ]  d^dTldC  (4.20) 


Niunerical  Integration 

Equation  4.20  must  be  integrated  numerically  by  a  Gauss 
Quadrature  formula  as  shown: 
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I  =  f(k.  T],  C)  didr\d(i 


"  E E E  ^j'  Cjt) 

i=l  j=l  Jc=l 

n  is  a  specified  integer,  and  from  this  choice  the  w's  and 
(5,/  »?j^  set.  n  =  2  integrates  a  cubic  equation 

exactly,  and  this  typically  gives  enough  accuracy.  However, 
n  =  4  is  used  in  the  stiffness  matrix  computations  because 
of  the  special  needs  of  the  unconfined  flow  algorithm  dis¬ 
cussed  later. 

Decomposition 

The  9  by  9  stiffness  matrix  is  next  reduced  to  an 
8  by  8  matrix  by  eliminating  terms  related  to  the  ninth  or 
internal  node.  This  is  done  by  partitioning  the  9  by  9 
equation , 

[iC]9{i2}9  =  {Q}g 

as  follows: 

{K}1  )[  ^9  J  [  0 

Here  Q,  =  0  because  no  flow  is  considered  to  enter  at  the 
internal  node.  This  matrix  equation  represents  two  equa¬ 
tions  where  the  second  one  can  be  solved  for  h,  and  substi¬ 
tuted  into  the  first  one.  This  gives 
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and 


h,  =  --^{K}l{h}, 

Agg 

So  the  actual  stiffness  matrix  [K]^  used  for  further  compu¬ 
tations  is 

[K]^  =  [A]g  -  -^{K}^{K}1  (4.21) 

Agg 

Discharge  Velocity 

The  discharge  velocity  {v}  for  each  element  can  be 
computed  by  the  general  version  of  Equation  4.8  as  follows: 

{v}  =  -[k]  [B]  (4.22) 

For  the  element  type  used  in  this  application,  discharge 
velocity  is  typically  computed  at  the  (i ,  ri ,  C)  =  (0,  0,  0) 
point  which  is  close  to  the  centroid  of  the  element.  The 
element  discharge  velocity  {v)^  at  this  point  thus  becomes 

(Wo  =  -Ik]  [B]Jh}^  (4.23) 

where  [B]q  is  the  [B]  matrix  evaluated  at  ($,  r?,  C)  =  (0,  0, 
0)  . 

Assemblage  with  Boundary  Conditions 
With  the  element  stiffness  matrices  computed,  one  can 
now  assemble  a  global  set  of  equations 
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=  {Q}a 


(4.24) 


{h}j,  is  the  set  of  global  heads  at  the  active  nodes  that 
must  be  computed.  If  the  boundary  condition  of  specified 
head  is  applied  to  a  node,  it  is  not  active  because  head  is 
already  known.  [K]^  is  computed  by  assembling  the  element 
stiffness  matrices  taking  into  account  the  active  nodes  as 
follows : 

[K]^  =  E 

J=1 

where  is  the  number  of  elements  and  [Kj^.  is  the  ith  ele¬ 
ment  stiffness  matrix  as  given  in  Equation  4.21. 

{Q}|.  is  the  combination  of  specified  flows  at  nodes 
plus  those  computed  from  a  specified  normal  discharge  velo¬ 
city  v^.  Let  {Q)^  be  flows  computed  at  four  nodes  of  a  face 
of  an  element  in  (^,  r))  computational  space  where  v^^  is 
applied.  Also,  let  the  shape  functions  {N}^  now  be  given  by 

(1  +  i  =  1,  2,  3,  4 

Then 


{C>}4  =  f  f  v^/eG  -~F^  {N},  d^dTi  (4.25) 
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where 


E  = 

F  = 

G  = 


a{r}^  a{r} 

di  di 

d{r}'^d{r} 

di  an 

a{r}^a{r} 

an  an 


Now  (r)  depends  only  on  four  nodes  on  th'  surface.  Let 
these  (X,  y,  z)  coordinates  be  designated  by  (x)^,  {y}^,  and 


where 


respectively. 

Then  {r} 

is 

simply 

{r}  = 

[i^]4[r]4 

'imi 

{0}J 

{0}!' 

[N],  = 

{Oil 

{ml 

{0)1 

[r]4  = 

{y}4 

,{0}! 

{Oil 

Also,  {0},^  is  a  zero  vector  of  length  4.  The  expressions 
for  E,  F,  and  G  are  now  easily  evaluated.  For  instance 


E  = 


[r] 


T-diml  d[N]^ 

"  a^  a$ 


All  the  terms  in  Equation  4.25  can  now  be  evaluated,  so  the 
expression  can  be  numerically  integrated  as  before.  Once 
the  {Q)^  are  computed,  they  are  assembled  into  the  global  Q 
vector  (Q)c- 

All  is  now  ready  to  solve  Equation  4.24  for  the  unknown 
heads  {h}^.  This  is  a  symmetric  banded  system  of  simulta¬ 
neous,  linear  equations  that  can  now  be  solved.  Once 
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the  heads  are  known  for  each  node  the  discharge  velocities 
can  be  computed  for  each  element  using  Equation  4.22.  If 
the  problem  is  a  confined  flow  problem,  the  solution  is  now 
complete.  However,  for  an  unconfined  flow  problem  this 
represents  only  the  initial  iteration.  Details  regarding 
unconfined  flow  problems  are  given  next. 

Unconfined  Flow  Problem 

The  unconfined  flow  problem  requires  an  iterative  pro¬ 
cess  to  determine  a  solution  because  of  the  free  surface 
being  unknown.  We  will  now  examine  this  problem  in  detail. 

Description  of  the  Problem 

To  see  the  complexity  in  two  dimensions,  consider  the 
earth  dam  shown  in  Figure  28.  Line  segment  DG  is  the  free 
surface  that  is  unknown.  The  exit  point  is  point  G,  and 
segment  GH  is  the  surface  of  seepage.  Notice  that  point  G 
does  not  necessarily  coincide  with  point  H  where  the  tail 
water  starts.  Line  segment  BCD  is  an  equipotential  line. 


E  F 
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and  since  DG  is  a  flow  line  at  steady-state  conditions,  DG 
is  perpendicular  to  BCD  at  D.  HIJ  is  an  equipotential  line, 
and  DEFG  is  removed  from  the  problem.  Points  A  and  B  are 
moved  far  enough  left  to  where  the  effects  of  the  dam  are 
not  felt,  and  the  same  is  done  to  points  J  and  K,  except 
they  are  moved  far  to  the  right. 

The  exit  point  becomes  an  exit  line  in  three  dimen¬ 
sions.  Further,  in  an  aquifer  with  wells,  there  are  any 
number  of  independent  tailwater  levels,  exit  lines,  and 
surfaces  of  seepage.  Nevertheless,  it  is  the  goal  of  this 
dissertation  to  allow  the  user  to  solve  all  these  problems 
without  restriction  of  grid.  It  is  also  extremely  tedious 
to  have  to  give  an  initial  guess  of  the  free  surface,  espe¬ 
cially  when  there  could  be  any  number  of  them  in  3-D,  so  no 
initial  guess  is  to  be  required. 

Geometric  Considerations 

Figure  29  shows  a  portion  of  an  FEM  grid  with  the  free 
surface  ABCDE  cutting  through  it.  Point  A  is  at  the  head¬ 
water  level,  and  point  F  is  at  the  tailwater  level. 

Reshaped  Elements 

First,  the  free  surface  cuts  the  grid  so  that  some 
elements  are  completely  out  of  the  problem  after  a  few  iter¬ 
ations,  and  others  have  new  shapes.  Some  become  very  small 
or  skewed,  while  others  have  five  sides.  In  three  dimen¬ 
sions  these  different  cases  become  even  more  complicated. 
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A 


Figure 

Clearly,  an  adaptive  grid  scheme  or  another  procedure  is 
more  advantageous . 

Exit  Point 

The  exit  point  E  also  causes  a  problem.  Points  B,  C, 
and  D  can  be  considered  to  have  no  flow  entering  at  these 
points  since  the  free  surface  is  a  flow  line.  Previous  work 
has  also  given  this  no-flow  condition  to  the  exit  point  E. 
This,  however,  is  incorrect,  and  if  done  this  way,  leads  to 
a  wrong  exit  point.  To  understand,  consider  Figure  30  where 
the  arrows  show  the  direction  of  flow.  Flow  is  lumped  at  a 
node  by  attributing  to  that  node  the  flow  in  the  vicinity  of 
that  node.  Although  no  flow  is  added  to  point  E  from  the 
top,  there  is  definitely  some  flow  from  the  surface  of  seep¬ 
age  part  below  the  node.  This  situation  becomes  even  more 
complex  in  three  dimensions.  Clearly,  it  would  be  better  to 


E 


60 


E 


Figure  30.  Exit  Point  E 

find  an  algorithm  that  completely  avoids  these  problems,  at 
least  during  the  iteration  process. 

Computational  Procedure 

After  trying  several  algorithms,  the  one  best  fulfill¬ 
ing  the  goals  of  this  dissertation  is  now  outlined: 

1.  Compute  and  store  for  each  element. 

2.  Assemble  [K]^  and  [KJ^q. 

3.  Solve  [K]go{h}j.o  =  {Q}^^  for 

4.  Compute  {Q}^  =  [K]^{h)j.o. 

For  j  =  1 ,  2 ,  3 ,  ... 


5. 

Check  for  convergence. 

6. 

Compute 

[K]gj  for  each  element. 

7  . 

Compute 

=  [K]^j.{h)^_..,  and 

assemble 

8. 

Compute 

different  global  stiffness  matrix 

9. 

Solve  [K],..{Ah},^  =  -{AQ)g.  for 

10. 

Update 

Ihicj  =  Ihlj,,.,  +  Uh),j. 
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End  j 

11.  Compute  the  final  position  of  the  free  surface. 

12.  Compute  the  final  element  discharge  velocities. 
Each  part  will  now  be  described  in  detail. 

Compute  and  Store  [K]^  for  Bach  Element 

The  first  thing  that  is  done  is  to  compute  the  stiff¬ 
ness  matrices  for  each  element.  Although  possibly  modified 
in  subsequent  iterations,  the  initial  material  properties 
(permeabilities)  are  used  here.  The  element  stiffness 
matrices  are  stored  and  used  later  in  the  iteration  process 

Assemble  [K]^  and  [K]^ 

The  unmodified  stiffness  matrix  [K]^  is  assembled  by 
using  both  active  and  inactive  nodes  in  the  assemblage  pro¬ 
cess.  Another  way  of  saying  this  is  that  the  boundary  con¬ 
ditions  are  not  considered.  is  determined  by  either 

assembling  the  global  stiffness  matrix  using  only  active 
nodes  (those  where  heads  are  not  specified)  or  modifying 
[K]^j  for  boundary  conditions. 

Solve  [K]go(*'>go  =  (Q>co 

This  is  a  standard  system  of  banded,  simultaneous, 
linear  equations  that  are  solved  for  the  heads  at  each  node 
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compute  {Q}^  =  [K]„{h>g<, 

The  flows  at  internal  nodes  are  typically  zero.  One 
exception  is  when  a  source  or  sink  is  specified.  The  nodes 
where  flow  is  specified  will  give  the  same  answer.  The 
nodes  where  head  is  specified  will  have  their  unknown  flows 
computed  by  this  operation. 

Check  for  Convergence 

The  basic  idea  is  to  iterate  until  the  free  surface  has 
stopped  moving.  Various  criteria  can  be  used  to  accomplish 
this  with  the  one  that  proved  best  being  now  presented. 
First,  compute  a  small  number  €, 

e  =  0.001  (z^  - 

where  z^^  is  the  maximum  z,  and  z^^^  is  the  minimum  z.  Next, 
compute  the  maximum  change  in  head  for  this  iteration  for 
all  nodes  that  are  still  in  the  flow  region  as  follows: 

A/Ux  =  max  ,)  J 

i 

where  i  is  the  ith  eligible  node.  When 

for  two  straight  iterations,  the  solution  has  converged. 

The  reason  two  iterations  are  used  is  that  for  some  rare 
problems  the  criteria  must  be  applied  twice  in  a  row  to 
achieve  true  convergence.  Otherwise,  false  convergence 
occurs  after  two  or  three  iterations. 
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Compute  [K]gj  for  Each  Element 

As  shown  in  Figure  29,  the  free  surface  is  above  some 
elements,  goes  through  others,  and  is  completely  below  yet 
others.  For  elements  that  are  completely  inside  the  flow 
region,  [K]gj  is  read  from  the  disc  and  used  as  is.  For 
those  completely  above  the  region,  the  strategy  of  using  a 
reduced  element  stiffness  matrix  is  used.  This  can  be 
either  in  steps  or  all  at  once.  The  reduction  factor  of 
0.001  is  used  here.  This  is  equivalent  approximately  to 
saying  that  no  flow  is  in  the  element.  Figure  31  shows  a 
2-D  example  of  what  is  done  for  elements  where  the  free 


Figure  31.  Free  Surface  through  2-D  Element 

surface  crosses  the  element  when  two  interpolation  points  in 
each  direction  are  used.  Basically,  the  numerical  integra¬ 
tion  process  in  this  case  involves  integrating  over  a  dis¬ 
continuity.  The  pressure  head  is  computed  at  each  integra¬ 
tion  point  used  in  integrating  Equation  4.20,  and,  if  less 
than  zero,  the  reduced  value  of  permeability  tensor  is  used. 
Also,  a  higher  number  of  integration  points  are  used  since  a 
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discontinuity  is  involved.  So  the  3-D  implementation  has 
the  total  head  computed  at  each  of  4  by  4  by  4  numerical 
integration  points  as 

h{^^,  Cj,)  =  Ci,))^{h}, 

where  i,  j,  and  k  range  from  1  to  4 .  The  pressure  head  hp 
for  an  arbitrarily  specified  datum  h^  is 


compute  {AQ}^-  and  {AQ}^. 

For  an  internal  node  the  flow  should  be  zero.  If  it  is 
not,  this  indicates  that  flux  is  going  across  the  boundary 
at  that  node,  and  it  is  part  of  a  moving  free  surface.  The 
technique  used  is  to  first  compute  element  flows  {AQj^j  by 

and  then  assemble  them  into  a  change  in  flow  for  each  node 
designated  by  {AQ}j,j.  For  an  external  node  where  the  flow 
is  specified,  the  change  in  flow  AQ'  becomes 

AQ'  ^  AO  -  Qe 

Here,  AQ  is  the  value  from  the  array  (AQ)(,j  for  that  node. 

compute  Different  Global  Stiffness  Matrix  [K]Q,j 

The  changes  in  head  {Ah)^^  could  be  computed  using  the 
equation 
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{A. 26) 


where  [K]gj  is  the  global  stiffness  matrix  assembled  from 
the  [K]gj.  However,  since  a  steady-state  solution  is 
desired,  it  is  expedient  to  employ  the  aerospace  engineers' 
strategy  that  states  for  a  numerical  expression  solving  for 
Ah 

(NUMERICS)  (Ah)  =  PHYSICS 

Now  the  physics  require  that  the  right-hand  side  of 
Equation  4.26  go  to  zero.  But  the  numerics  of  the  left-hand 
side  can  be  simplified  to  achieve  this  goal  by  using  the 
unmodified  stiffness  matrix  [K]^  (which  has  already  been 
computed)  and  applying  the  current  boundary  conditions  to 
obtain  a  different  global  stiffness  matrix  [K]g,j.  As  will 
now  be  explained  using  Figure  32,  the  only  place  where 


Figure  32.  Surface  of  Seepage 
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boundary  conditions  change  is  on  the  surface  of  seepage. 

Line  segment  AB  is  an  exit  line  on  an  exit  face,  which  is 
shown  only  with  the  external  sides  of  the  brick  elements  at 
the  exit  face.  Each  node  below  the  line  of  seepage  and 
above  or  on  the  tailwater  area  will  have  a  negative  flow, 
since  flow  is  leaving  the  system  at  that  point.  So  when 
such  a  node  is  encountered,  its  pressure  head  is  zero,  and 
its  boundary  condition  for  head  is 

h  =  z  - 

Now,  if  in  the  course  of  iteration  the  flow  for  this  node 
becomes  positive,  then  the  free  surface  has  fallen  below 
this  node,  and  its  boundary  condition  must  be  changed.  What 
is  done  is  illustrated  in  the  2-D  sketch  shown  in  Figure  33 
where  the  region  above  the  exit  point  is  made  impervious. 

So  the  node  with  the  status  change  will  have  its  boundary 
condition  shifted  from  specified  head  to  the  no-flow  bound¬ 
ary  condition. 


Figure  33.  2-D  Exit  Point  Analysis 
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A  node  can  also  be  in  the  no-flow,  impervious  condition 
and  need  to  change  back  to  a  node  below  the  exit  line.  This 
is,  in  fact,  a  case  of  a  rising  exit  line.  What  is  done  is 
to  check  the  currently  impervious  node  on  the  surface  of 
seepage  for  the  condition 

h  z  -  hjy 

If  this  condition  is  passed,  then  the  node  is  changed  back 
to  a  specified  head  node  with  the  boundary  condition 

h  =  z  - 

again  applied. 

Compute  (Ah}gj 

Let  the  unmodified  stiffness  matrix  which  has  now 

been  modified  at  iteration  j  by  the  current  boundary  condi¬ 
tions  be  designated  as  Then  the  actual  system  of 

equations  used  to  compute  the  change  in  head  at  iteration  j 
is 

(4.27) 

Equation  4.27  is  a  symmetric,  positive  definite,  banded 
system  of  simultaneous  linear  equations  that  can  now  be 
solved  for  {Ahj^j. 

Update  {h}gj 

The  active  nodes  can  now  have  their  total  heads  updated 
as  follows: 
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{h}^^  =  {iJlc.j.i  +  (4.28) 

One  iteration  is  now  complete,  and  the  cycle  is  repeated. 

The  process  is  continued  until  convergence  is  achieved. 

Determine  Final  Position  of  Free  Surface 

The  determination  of  the  final  position  of  the  free 
surface  for  a  structured  grid  algorithm  is  rather  straight¬ 
forward  .  However,  the  thrust  of  this  work  is  to  do  a  free 
surface  problem  without  any  required  structure  at  all  to  the 
grid.  The  determination  of  the  free  surface  in  this  case  is 
significantly  more  difficult,  especially  at  the  exit  line. 
First,  consult  Figure  34  to  understand  the  complexities  away 
from  the  exit  line.  Here  two  tetrahedral  elements,  ABEC  and 
AGED,  of  an  unstructured  3-D  grid  are  shown  with  the  free 
surface  denoted  by  the  cross-hatched  plane  crossing  these 
elements.  First,  the  reason  the  free  surface  crosses  rhece 
elements  is  that  the  pressure  head  at  point  A  is  less  than 
zero,  and  at  points  B,  C,  D,  and  E  the  pressure  head  is 
greater  than  zero.  Now  it  is  best  to  keep  the  same  data 
structure  (FEM  mesh)  for  the  purpose  of  computing  discharge 
velocities  at  the  elements  and  providing  output  to  post¬ 
processor  programs.  So  one  strategy  is  to  "clip"  an 
affected  element  by  the  free  surface  to  create  two  elements. 
Then  one  element  would  be  completely  above  the  free  surface, 
and  one  would  be  completely  below  the  free  surface.  How¬ 
ever,  this  leads  to  significant  complications.  Notice, 
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Figure  34.  Free  Surface  across  Grid 


first  of  all,  that  no  problem  exists  in  Figure  34  where  each 
four-sided  piece  is  changed  into  a  four-sided  piece  and  a 
five-sided  piece.  However,  consider  the  effect  of  a  six- 
sided  element  being  clipped  by  the  free  surface  with  the 
piece  left  below  the  free  surface  shown  in  Figure  35.  What 
is  left  is  a  seven-sided  piece,  no  longer  conforming  to  the 
allowable  shapes  of  the  eight-node  brick  element.  This 
problem  could  be  alleviated  by  dividing  the  seven-sided 
figure  into  acceptable  pieces.  However,  the  variety  of 
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Figure  35.  Seven-Sided  Piece 


unacceptable  shapes  warranted  the  search  for  another  solu¬ 
tion.  Point  A  in  Figure  34  with  a  negative  pressure  head 
was  moved  to  a  point  on  the  free  surface  where  the  pressure 
head  is  zero.  This,  however,  creates  a  new  dilemma.  Any 
point  in  the  cross-hatched  surface  representing  the  free 
surface  has  zero  pressure  head,  so  some  criteria  must  be 
devised  to  decide  which  point  will  be  chosen. 

After  trying  different  algorithms,  the  one  that  works 
well  and  was  implemented  will  now  be  given: 

I.  Identify  a  node  (such  as  point  A  in  Figure  36)  with 
pressure  head  less  then  zero. 

A.  Find  a  line  segment  emanating  from  A  so  that  the 
second  point  of  the  line  segment  (such  as  point  B 
in  Figure  36)  has  pressure  head  greater  than  zero. 

1.  Compute  the  value  of  the  parameter  s 


0  i  s  i  1 


with  s  =  0  at  point  A  and  s  =  1  at  point  B 
where  the  pressure  head  is  zero  (point  O  in 
Figure  36) .  Linear  interpolation  gives 


A 


B 

Figure  36.  Line  Segment  AB 


•®o 


(4.28) 


where  h  .  is  the  pressure  head  at  point  A, 
hpg  is  the  pressure  head  at  point  B,  and  s^ 
is^the  value  of  s  at  point  O. 

2.  Compute  the  new  position  of  point  A  along  AB  so 
it  resides  on  the  free  surface  by 

{r}^  =  (1  -  s^)  {r}^  +  sjr}^ 

where  {r)^  is  the  new  position  of  point  A,  {r}^ 
is  the  original  position  of  point  A,  and  {r)g 
is  the  position  of  point  B. 

3.  Compute  the  distance  squared  from  point  A  to 
point  O.  Save  that  number  with  the  current 
position  of  the  tree  surface. 

B.  Repeat  step  A  for  other  line  segments  starting 

with  point  A.  Keep  the  one  (rl^  with  the  smallest 
distance  squared  value. 
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II.  Process  all  other  node  points  just  above  the  free 

surface  that  need  to  be  moved  in  the  same  way  as  in 
step  I . 

A  node  point  on  the  surface  of  seepage  cannot  be  done 
exactly  this  way.  The  reason  is  that  the  value  of  s^  can¬ 
not  be  computed.  To  understand,  first  consider  the  2-D  case 
illustrated  in  Figure  37.  Line  segment  AB  is  a  line  on  the 
surface  of  seepage,  and  E  is  the  exit  point  to  be 
computed.  Now  AB  contains  an  exit  point  because  (1)  the 
nodes  have  surface  of  seepage  boundary  conditions  and 
(2)  node  A  was  flagged  as  an  impervious  node  and  node  B  was 
flagged  as  a  specified  head  node  at  the  last  iteration  cycle 
before  convergence.  Points  a,  b,  and  c  are  the  three  most 
adjacent  free  surface  nodes  that  have  been  moved  from  their 
original  position  to  their  new  positions,  respectively,  to 
be  used  for  computing  point  E. 


A 


B 

Figure  37.  Exit  Point  Computation 

Exit  point  E  must  be  determined  by  extrapolation.  What 
is  used  is  a  difference  formulation  using  points  1  through 
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3.  First,  the  slope  at  point  a  is  computed  by  the  backward 
difference  expression 


(Vb  -  ya> 

(Xt,  -  xj 


In  a  similar  manner  the  slope  at  point  b  is 


^b 


(Vc  -  Vb^ 

(x^  -  x^) 


From  these  results  the  change  in  slope  c^  can  be  computed  by 


c 


m 


(^b  - 

i^b  - 


We  can  now  compute  the  slope  at  point  c  using 

With  x^,  y^,  and  m^  all  known  this  line  can  be  intersected 
with  line  segment  AB  to  obtain  the  intersection  point  E. 

This  is  accomplished  by  substituting  into  the  equation 

y  =  +  m^(x  -  x^) 

the  parametric  equation 

{r}  =  (1  -  s) {r}^  +  ^{r}^ 

and  solving  for  s.  If  s  does  not  lie  between  zero  and  one, 
the  intersection  point  occurs  outside  the  range  of  line 
segment  AB,  and  point  B  is  taken  as  the  exit  point. 

The  3-D  case  is  much  more  difficult.  Figure  38  shows  a 
plan  view  of  a  portion  of  the  free  surface  just  computed 
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with  the  exit  line  FG  needing  to  be  computed.  What  is  sig¬ 
nificantly  more  difficult  is  the  need  to  extrapolate  off  an 
unstructured  surface  to  get  the  intersection  of  that  surface 
with  the  surface  of  seepage  yielding  the  points  on  the  exit 
line  FG.  After  some  analysis  of  these  difficulties,  it  was 
decided  to  take  another  approach.  Figure  39  shows  an  almost 
vertical  exit  face  where  it  has  been  determined  that  point  A 
must  move  down  on  this  surface  to  become  part  of  the  exit 
line.  Now  at  the  last  iteration,  point  B  has  the  surface  of 
seepage  boundary  condition 

h  =  z  - 

and  point  A  is  denoted  as  impervious  with  a  computed  head 
value 

z  ~  ho 
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A 


B 


Figure  39.  Exit  Face 


The  amount  point  A  must  move  down  line  segment  AB  is  deter¬ 
mined  by  computing  a  value  for  the  parameter  s.  First, 
compute  Sg  by 


So  = 


Sfi  ~  ^A 


(4.29) 


The  final  value  of  s  is  computed  from 

=  (1  -  y)  Sg  +  Y  (4.30) 

with 

0 . 5  Y  ^  1 

From  experience  the  value  used  in  this  work  is  0.8.  The  new 
position  of  point  A  on  AB  is  computed  as  before 

{r}  1  =  (1  -  Sj)  {r}^  +  {z}  g 

As  in  the  other  free  surface  computations,  point  A  can 

move  in  more  ways  than  along  line  segment  AB.  In  Figure  39, 
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for  instance,  line  segments  AC  and  AD  must  also  be  consid¬ 
ered.  As  before,  the  line  segment  is  kept  with  the  smallest 
distance  of  movement  of  point  A  from  its  original  position. 

Compute  Final  Element  Discharge  Velocities 

Now  that  the  nodes  near  the  free  surface  have  been 
moved  to  be  on  the  free  surface,  the  element  data  must  be 
modified  so  that  all  elements  are  either  completely  above 
the  free  surface  or  completely  below  it.  Figure  40  shows  an 
originally  square  mesh  with  the  new  shapes  of  the  elements 
conforming  to  the  free  surface  shown.  Those  elements  above 
the  free  surface  will  be  given  a  value  for  discharge 
velocity  of  zero,  and  those  below  the  free  surface  will  have 
their  discharge  velocities  computed  using  Equation  4.23. 
Additional  care  must  also  be  taken  to  ensure  that  the  newly 
formed  elements  that  become  skewed  are  handled  properly. 

They  can  be  either  removed  fi'om  the  grid  or  kept  and  given 
zero  discharge  velocity.  In  this  work  they  are  kept  and 
given  zero  values  of  discharge  velocity  for  compatibility 
with  postprocessing  programs. 

Achieving  the  newly  formed  grid  as  shown  in  Figure  40 
is  more  difficult  than  it  first  appears.  First,  consider 
the  2-D  case  where  Figure  41  shows  the  three  steps  used  in 
the  process.  Step  1  consists  of  determining  the  position  of 
the  free  surface  nodes  on  the  proper  line  segments  (the  dots 
in  Figure  41).  In  step  2  the  free  surface  nodes  are  moved 
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to  their  new  positions.  Step  3  is  the  most  complicated  in 
that  additional  movement  of  nodes  is  required  to  get  ele¬ 
ments  either  completely  below  or  completely  above  the  free 
surface . 

The  3-D  version  of  this  problem  is  even  more  compli¬ 
cated.  This  is  because  of  the  many  different  shapes  the  3-D 
brick  element  can  take  as  shown  in  Figure  27.  Also,  because 
of  the  general  nature  of  the  3-D  element,  establishing  what 
is  on  top  and  what  is  on  bottom  becomes  blurred.  However, 
as  illustrated  in  Figure  42,  a  3-D  version  of  the  algorithm 
presented  above  was  developed  as  part  of  this  work.  Fig¬ 
ure  42  shows  a  brick  element  with  free  surface  nodes 


3  2 


3 


2 


Before 

Figure  42. 


After 

3-D  Version 
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occur ing  exactly  at  nodes  3  and  8.  The  first  plot  shows  the 
element  before  any  change,  and  the  second  plot  shows  the 
element  after  it  was  modified  to  fit  the  free  surface.  Step 
3  of  the  algorithm  collapsed  node  7  down  to  node  8.  In 
fact,  the  first  step  3  algorithm  considered  checks  each  of 
the  six  faces  of  an  element  to  see  if  any  two  nodes  on  a 
diagonal  of  a  face  are  bath  free  surface  nodes.  If  two 
exist,  then  any  other  node  on  the  face  having  a  head  value 

h  <  z- 

will  be  collapsed  to  one  of  the  diagonal  nodes. 

This  algorithm,  however,  does  not  always  work.  The  2-D 
version  of  this  algorithm,  for  instance,  works  on  all  verti¬ 
cal  columns  of  six  elements  each  in  Figure  40  except  the 
last  one.  The  isolated  part  of  the  grid  where  the  problem 
exists  is  shown  in  Figure  43.  Figure  44  shows  the  three 
step  sequence  and  how  it  breaks  down.  The  middle  element 
marked  by  the  X  does  not  have  two  free  surface  nodes 
on  a  diagonal,  so  the  third  step  is  not  completed.  There¬ 
fore,  a  higher  order  iterative  algorithm  must  be 
implemented . 

The  following  algorithm  was  developed  as  part  of  this 

work ; 
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Repeat  until  no  change  occurs 

For  n  =  1  to  the  number  of  elements 

If  the  free  surface  intersects  the  element,  then 

For  i  =  1  to  12  /*  Test  all  line  segments.  */ 

If  a  free  surface  node  is  connected  to  a  node 
above  the  free  surface,  then 

Change  the  coordinates  and  head  of  the 
node  above  the  free  surface  to  that  of  the 
free  surface  node. 

Flag  the  moved  node  as  a  free  surface 
node,  but  as  an  artificially  generated 
one . 

End  If 

End  i 

End  If 

End  n 

End  Repeat 

The  algorithm  can  be  utilized  in  two  ways: 

1.  Adjust  only  those  nodes  immediately  above  the  free 
surface  and  leave  the  other  nodes  alone. 

2.  Adjust  those  nodes  immediately  above  the  free 
surface  and  have  the  nodes  farther  above  the  free 
surface  to  collapse  to  the  free  surface. 

Appendix  A  contains  a  FORTRAN  listing  and  description  of  the 

free  surface  algorithms  described  in  this  chapter  with 

implementation  of  Option  2  for  utilizing  the  algorithm. 
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CHAPTER  V 


THREE-DIMENSIONAL  SEEPAGE  PACKAGE  AND  GROUNDWATER  MODEL 

Introduction 

A  three-dimensional  seepage  package  and  groundwater 
model  was  developed  that  concentrated  on  three  areas: 

1.  Sophisticated  numerics  to  allow  both  structured  and 
unstructured  grids  with  a  minimum  of  input  for  free 
surface  problems. 

2.  Interface  to  state-of-the-art  grid  generation 
technology. 

3.  Interface  to  state-of-the-art  scientific 
visualization  technology. 

The  numerics  were  discussed  in  detail  in  Chapter  IV.  This 
chapter  describes  the  components  of  the  3-D  model. 

Grid  Generation 

EAGLE,  written  originally  for  aerospace  engineering 
applications  (computational  fluid  dynamics  (CFD))  to  gener¬ 
ate  structured  grids  for  finite  volume  solvers,  was  used  as 
the  primary  grid  generation  tool  in  this  work.  It  has 
extensive  algebraic  grid  generation  capabilities,  as  well  as 
state-of-the-art  elliptic  smoothing  capabilities.  However, 
since  the  FEM  was  the  computational  tool  used,  partly 
unstructured  of  totally  unstructured  grids  may  also  be  used 
with  the  developed  computer  program. 


Conversion  to  FEM  Format 

The  finite  volume  grid  from  EAGLE  with  potentially 
several  blocks  must  be  converted  to  the  final  data  needed  by 
the  FE  program.  Five  things  are  involved: 

1.  Getting  title,  datum,  and  material  property 
information. 

2.  Applying  boundary  conditions. 

3.  Combining  the  different  blocks  into  one  FEM  grid. 

4.  Applying  bandwidth  minimization. 

5.  Writing  an  output  file  containing  this  information. 
Appendix  B  contains  a  description  of  the  program  written  to 
accomplish  these  tasks. 


Getting  Data 

Getting  title,  datum,  and  material  property  information 
is  a  simple  matter  of  prompting  the  user  for  this 
information. 


Applying  Boundary  Conditions 
The  EAGLE  finite  volume  grid  being  structured  makes 
applying  boundary  conditions  much  easier  than  usual.  What 
is  done  is  to  ask  for  the  following  data: 

I  BLOCK.  II,  Jl,  Kl.  12,  J2,  K2 ,  IBC,  BV 
where  the  meanings  of  the  variables  are  as  follows: 
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I BLOCK  -  Block  number. 


11  -  First  i  value. 

J1  -  First  j  value. 

K1  -  First  k  value. 

12  -  Second  i  value. 

J2  -  Second  j  value. 

K2  -  Second  k  value. 

IBC  -  Boundary  value  code. 

BV  -  Boundary  value. 

The  possible  values  of  IBC  and  BV  are  given: 


IBC  Value 


1 

2 

12 


-1 


Boundary  Value 

Specified  head. 

Exit  face.  Head  is  computed 
hY  h  =  z  -  . 

Border  of  exit  face  and  tail- 
water  level.  Head  is 
computed  the  same  as  IBC  =  2. 

Specified  flow. 


-2 


Specified  discharge  velocity. 


What  this  allows  the  user  to  do  is  to  apply  a  boundary 
condition  on  an  entire  face.  For  example,  consider  the 
guadrilateral  earth  dam  represented  by  one  block  shown  in 
Figure  45.  The  i  index  of  the  block  varies  in  the  horizon¬ 
tal  direction,  the  j  index  varies  into  the  paper,  and  the  k 
index  varies  along  lines  such  as  ABE  and  IG.  Suppose 
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Figure  45.  Quadrilateral  Earth  Dam 


this  is  a  101  by  101  by  101  points  block.  Also,  let  the  (x, 
y,  z)  coordinates  and  (i,  j,  k)  values  be  as  follows: 


Letter 

X 

Y 

Z 

i 

i 

k 

A 

0 

0 

0 

1 

1 

1 

B 

80 

0 

12 

1 

1 

81 

C 

0 

200 

0 

1 

101 

1 

D 

80 

200 

12 

1 

101 

81 

E 

100 

0 

15 

1 

1 

101 

F 

100 

200 

15 

1 

101 

101 

G 

110 

0 

15 

101 

1 

101 

H 

110 

200 

15 

101 

101 

101 

I 

210 

0 

0 

101 

1 

1 

J 

210 

200 

0 

101 

101 

1 

headwater 

value 

of  12 

ft  can 

now  be 

specified  by 

data : 
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I BLOCK 

11 

J1 

K1 

12 

J2 

K2 

IBC 

BV 

1 

1 

1 

1 

1 

101 

81 

1 

12 

In  a  similar 

manner. 

the  boundary 

conditions  of 

the  exit 

face  are 

defined  by 

the  two 

lines 

of  data: 

I BLOCK 

11 

Jl 

K1 

12 

J2 

K2 

IBC 

BV 

1 

101 

1 

2 

101 

101 

101 

2 

— 

1 

101 

1 

1 

101 

101 

1 

12 

- 

The  -  for  Bv  indicates  that  no  data  is  needed,  but  a  zero 
(0)  must  be  supplied  for  the  sake  of  completeness  of  data. 


Combining  Different  Blocks 

The  process  of  combining  different  blocks  into  a  single 
finite  element  mesh  consists  of  the  following: 

1.  Obtaining  a  single  set  of  consecutively  numbered 
nodes  and  elements. 

2.  Removing  duplicate  nodes  and  modifying  the  node 
numbering  and  element  connectivity  matrix  to 
reflect  the  changes. 

Note  that  duplicate  nodes  can  occur  in  an  O  type  grid  along 
the  cut  as  well. 


Applying  Bandwidth  Minimization 
There  are  different  optimization  schemes  with  varying 
degrees  of  sophistication.  Some  rearrange  the  node  numbers 
while  others  modify  the  element  numbering.  It  was  decided 
to  keep  the  element  numbers  the  same  in  this  program  to 
allow  the  ability  to  output  data  to  plotting  programs  like 
FAST.  Therefore,  a  bandwidth  minimization  algorithm  was 
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implemented  to  modify  the  node  numbers.  Because  memory  is 
plentiful  on  the  Cray,  a  straightforward  algorithm  was  used 
as  follows: 

1.  Create  an  adjacency  table  stating  what  other  nodes 
are  connected  to  a  given  node  and  how  many. 

2.  Place  the  new  node  1  at  various  places  on  the  grid. 
Then  number  the  adjacent  nodes  to  the  new  node  1  as 
2,  3,  4,  etc. 

3.  Go  to  the  new  node  2  and  number  the  adjacent  nodes 
to  that  node  which  have  not  been  numbered  as  5,  6, 
7,  etc.  Continue  this  process  until  all  the  nodes 
have  been  numbered. 

4.  After  trying  as  many  of  the  possible  positions  for 
the  new  node  1  as  desired,  keep  the  one  with  the 
smallest  bandwidth. 


Writing  an  Output  File 

The  last  thing  to  be  done  is  to  write  an  output  file 
containing  the  FEM  grid.  This  includes  all  the  preliminary 
data,  boundary  condition  information,  and  discharge  velocity 
data . 


Cray  Version  of  FEM  Program 
The  finite  element  program  was  converted  to  the  Cray 
YMP  with  as  much  voctorization  as  possible  being  accom¬ 
plished.  The  solution  of  the  banded  system  of  equations  was 
solved  by  using  an  out-of-core  solver  with  blocks  of  size 
1250.  The  forward  loop  of  the  Gauss  Elimination  algorithm 
was  vectorized,  but  the  back  substitution  loop  was  not. 
Nevertheless,  a  problem  having  10,644  nodes  that  took 
approximately  12  hr  on  the  Silicon  Graphics  Power  Series 


88 


4D/220  GT  took  3  min,  20  sec  on  the  Cray  YMP  6/128.  Even 
more  savings  can  be  attained  by  using  an  in-core  solver  and 
further  vectorization .  The  initial  savings  were  so  dramatic 
that  this  enhancement  will  be  done  at  a  future  date. 

Scientific  Visualization 

An  output  file  containing  the  results  was  written  to 
link  with  various  scientific  visualization  tools.  Of 
course,  this  file  differs  depending  on  which  tool  is  used. 
Appendix  C  contains  a  subroutine  that  writes  a  file  for  Flow 
Analysis  Solver  Toolkit  (FAST)  (Bancroft,  Kelaita,  McCabe, 
Merritt,  Plessel ,  Globus,  and  Semans  1991).  The  output  is 
similar  to  that  required  by  PLOT3D  or  Scientific  Visualiza¬ 
tion  Systems  (program  SSV)  sold  by  Sterling  Software.  How¬ 
ever,  the  data  must  be  converted  to  a  binary  file  written  in 
C  to  be  input  to  SSV  or  the  beta  version  of  FAST. 

Conversion  from  Element  Data  to  Node  Data 
The  discharge  velocities  are  computed  for  a  given  ele¬ 
ment  at  the  average  value  of  the  (x,  y,  z)  coordinates  of 
the  nodes  of  that  element.  However,  to  link  with  finite 
volume  scientific  visualization  tools,  values  are  needed  at 
each  node.  Consider  Figure  46  to  visualize  the  process  used 
to  get  node  values.  Discharge  velocity  at  node  n  is  com¬ 
puted  by  some  type  of  average  of  the  discharge  velocities 
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from  points  A,  B,  and  C.  The  program  in  Appendix  B  uses  a 
straight  average.  However,  another  popular  choice  is  a 
weighted  "one  over  distance  squared"  average  given  by 


where  d^^  is  the  distance  between  point  A  and  node  n,  dg^  is 
the  distance  between  point  B  and  node  n,  and  d;.^  is  the 
distance  between  point  C  and  node  n.  D  is  given  by 


Figure  46.  Element  to  Node  Conversion 


Conversion  from  Finite  Element  to  Finite  Volume  Format 
Now  that  values  of  head  and  velocity  exist  at  each 
node,  the  block  data  structure  must  be  reproduced  to  output 
a  results  file  for  FAST.  The  complication  is  that  node  num¬ 
bers  are  scrambled  from  bandwidth  minimization,  nodes  are 
eliminated  when  the  blocks  are  combined,  and  the  block  sizes 
are  discarded  when  creating  the  FEM  grid.  The  last  dilemma 
is  solved  by  writing  an  output  containing  the  block  sizes. 
Further,  since  the  element  numbers  have  been  preserved  in 
the  bandwidth  minimization  process  and  the  elements  were 
created  in  order  of  input,  the  element  connectivity  matrix 
can  be  used  to  write  the  output  file.  To  understand,  con¬ 
sider  Figure  47  showing  a  3-  by  3-  by  3-point  block.  First, 
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there  are  2x2x2  elements.  So  a  triple  loop  with  i2  =  j2 

=  k2  =  3  and  mstart  =  the  element  number  just  before  the 

block  begins  can  be  set  up  as  follows: 

mm  =  mstart 
For  k  =  1  to  k2-l 
For  j  =  1  to  j2-l 
For  i  =  1  to  i2-l 
mm  =  mm  +  1 

Output  results  for  the  1st  node  of  element  m. 
End  i 

/*  Do  last  one  on  i  row.  */ 

Output  results  for  the  2nd  node  of  element  m. 

End  j 

/*  Do  last  row  in  i,  j  plane.  */ 

mm  =  mm  -  i2  +  1 
For  i  =  1  to  i2-l 
mm  =  mm  +  1 

Output  results  for  the  4th  node  of  element  m. 

End  i 

Output  results  for  the  3rd  node  of  element  m. 

End  k 

/*  Do  last  i,  j  plane.  */ 

mm  =  mm  -  (i2-l)  *  (j2-l) 

For  j  =  1  to  j2-l 
For  i  =  1  to  i2-l 
mm  =  mm  +  1 

Output  results  for  the  5th  node  of  element  m. 

End  i 

Output  results  for  the  6th  node  of  element  m. 

End  j 

/*  Do  last  row  of  last  i,  j  plane.  */ 

mm  =  mm  -  i2  +  1 
For  i  1  to  i2-l 
mm  =  mm  +  1 

Output  results  for  the  8th  node  of  element  m. 

End  i 

Output  results  for  the  7th  node  of  element  m. 
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This  process  can  be  combined  into  a  single  set  of  loops 
by  use  of  a  variable  (IPICK  in  the  listing  below)  that 
selects  which  node  of  the  given  element  should  be  processed. 
A  FORTRAN  implementation  of  the  procedure  follows: 


C 

C 

C 

C 

C 

C 

C 


DIMENSION  IPICK(2,  2,  2) 

DATA  IPICK  /I,  2,  4,  3,  5,  6,  8,  7/ 

IROW  =12-1 

IPLANE  =  (J2  -  1)  *  IROW 

DO  K  =  1,  K2 
KK  =  1 

IF  (K.EQ.K2)  KK  =  2 

DO  J  =  1,  J2 
JJ  =  1 

IF  (J.EQ.J2)  JJ  =  2 

DO  I  =  1,  12 
II  =  1 

IF  (I.EQ.I2)  II  =  2 

MM  =  (K  -  KK)  *  IPLANE  +  (J  -  JJ)  *  IROW  +  I 
&  -  II  +  1  +  NSTART 

IPOS  =  IPICK(II,  JJ,  KK) 

Output  results  for  the  IPOS'th  node  of 
element  MM. 

END  DO 
END  DO 
END  DO 


The  subroutine  given  in  Appendix  C  uses  this  streamlined 
procedure . 

Visualization  Techniques 

It  turns  out  that  programs  such  as  FAST  written  for 
displaying  results  of  CFD  computations  using  structured 
grids  have  capabilities  very  useful  in  displaying 


seepage/groundwater  results.  One  of  the  goals  of  this  work 
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is  to  apply  aerospace  engineering  technology  to  the  flow  of 
groundwater  and  seepage  under  dams.  This  has  proven  to  be  a 
very  successful  endeavor.  Specific  tools  and  their  use  will 
now  be  described  with  actual  examples  given  in  Chapter  VI. 

Initial  and  Final  Grid 

Of  course,  the  first  thing  that  is  needed  is  the  dis¬ 
play  of  the  generated  grid.  Different  surfaces  of  the  grid 
are  selected  for  viewing  with  various  options  available, 
including  grid  lines,  continuous  tone  shading,  and  translu¬ 
cent  shading.  For  unconfined  flow  problems  the  grid  is 
modified  to  the  shape  of  the  free  surface,  and  the  part  of 
the  grid  where  there  is  no  longer  any  water  is  collapsed  to 
the  free  surface.  Thus,  it  is  extremely  helpful  to  view  the 
final  shape  of  the  grid,  as  well  as  the  original  grid,  to 
visualize  the  quality  of  the  solution. 

Isolevels 

One  important  aspect  of  the  flow  net  for  2-D  applica¬ 
tions,  such  as  the  one  given  in  Figure  12,  is  equipotential 
lines.  The  3-D  version  of  an  equipotential  line  is  a  sur¬ 
face  in  3-D  space  where  the  potential  is  a  constant  value. 

An  isolevel,  typically  used  to  look  at  where  a  particular 
value  of  pressure  exists  around  a  structure  such  as  an  air¬ 
foil,  is  simply  ideal  for  this  alternate  use.  Different 
isolevels  in  a  seepage  or  groundwater  problem  can  be  readily 
analyzed  by  the  practicing  engineer  to  check  the 
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reasonableness  of  the  solution.  Incorrect  boundary  condi¬ 
tions,  for  instance,  will  become  very  apparent. 

Color  Contours  on  a  Surface 

An  alternative  to  the  isolevels  is  color  contours  on  a 
given  surface.  The  surface  can  be  either  one  of  the  i-j , 
j-k,  or  k-i  planes  of  the  structured  grid  or  an  arbitrary 
cut  through  the  grid.  Successive  planes  can  also  be  viewed 
to  see  the  progression  of  results  as  well. 

Particle  Traces 

The  other  aspect  of  the  flow  net  given  in  Figure  13  is 
the  flow  lines.  While  it  is  not  possible  in  a  general  3-D 
problem  to  draw  a  flow  net,  the  particle  traces  typically 
used  by  the  aerospace  engineer  to  trace  air  flow  to  find, 
for  instance,  vortices  are  ideal  to  show  the  flow  of  water 
through  the  soil.  The  current  version  of  FAST  does  a 
steady-state  computation  of  a  particle  trace  which  works 
very  well  in  this  application. 
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CHAPTER  VI 


3-D  GROUNDWATER  FLOW  IN  AN  AQUIFER 

Introduction 

This  chapter  contains  a  three-dimensional  problem 
demonstrating  the  application  of  the  seepage/groundwater 
model  described  in  Chapter  V.  The  problem  that  is  given  is 
3-D  groundwater  flow  in  an  aquifer.  EAGLE  is  used  to  gener¬ 
ate  the  grids,  the  new  3-D  sc-^page, /groundwater  model  is  used 
to  do  the  computations,  and  FAST  is  used  to  display  the 
results.  Simpler  versia.::'  of  the  problem  are  first  done  to 
accomplish  the  following; 

1.  Compare  with  known  results. 

2.  Investigate  the  quality  of  the  selected  grid. 

3.  Illustrate  the  scientific  visualization 
capabilities  of  FAST. 

Of  specific  importance,  also,  is  whether  the  numerics  of  the 
unconfined  portion  of  the  implemented  algorithms  work 
properly . 


Groundwater  Flow  in  an  Aquifer 
The  problem  of  groundwater  flow  in  an  aquifer  will  now 
be  presented. 
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General  Description  of  Problem 
The  problem  of  groundwater  flow  in  an  aquifer  involves 
flow  of  water  subject  to  pumping  and  recharge.  Pumping 
occurs  through  wells  with  varyiri^  yields.  Natural  recharge 
occurs  by  filtration  through  the  bed  of  rivers  crossing  the 
aquifer  area  and  by  seepage  through  zones  of  limited  extent 
located  at  the  boundary  of  the  aquifer.  Artificial  recharge 
occurs  by  filtration  from  the  ground  surface.  Finally,  the 
aquifer  is  characterized  by  zones  with  significantly  differ¬ 
ent  material  types  with  some  having,  for  instance,  a  strong 
anisotrophy  of 


Simplified  Problem 

The  simplified  problem  will  now  be  discussed  and 
illustrated. 

Description  of  Simplified  Problem 

The  simplified  problem  is  shown  in  Figure  48  and  con¬ 
sists  of  a  very  small  homogeneous  aquifer  1,000  ft  long, 

500  ft  wide,  and  120  ft  deep.  Two  partially  penetrating 
wells  with  the  following  characteristics  have  been  placed  in 
the  aquifer: 

Well _ X _ Y _ Radius  Penetration _ Q _ 

1  250'  200'  1'  40'  100  cfm 

2  600'  250'  2’  80'  200  cfm 
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RIVER 


a.  Plan  View 


b.  Front  View 

Figure  48.  Aquifer  with  Two  Wells 
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Only  confined  flow  exists  in  the  aquifer,  and  it  is 
impervious  on  three  sides,  as  well  as  the  top  and  bottom. 

The  remaining  side,  however,  is  recharged  by  a  river  modeled 
by  a  constant  head  of  50  ft  above  the  top  of  the  aquifer. 

The  wells  are  modeled  by  the  flows,  and  Qj,  being 
distributed  uniformly  on  the  sides  of  the  respective  wells. 
The  permeability  of  the  soil  is  0.1  ft/min. 


Analytic  Solution  for  a  Partially  Penetrating  Well 

An  analytic  solution  to  a  partially  penetrating  well 
(Figure  49)  of  thickness  t^^,  penetration  b,  well  radius  r^, 
permeability  k,  and  flow  Q  has  been  developed  (Muskat  1946) , 
and  it  will  now  be  given.  First,  define  the  variables 


a  = 


2t. 


c  =  - 


2t. 


P  = 


oc. 


2  t. 


(6.1) 


where  (r,  z)  =  (0,  0)  is  located  at  the  top  of  the  aquifer 
and  at  the  center  of  the  well.  Next,  define  the  function 


Z(s,  y) 


oa 


n=0 


1 

(n  +  y) ^ 


(6.2) 


For  small  values  of  a. 
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^  ^r{i  -  c  -  (i)r(i  +  c  -  P) 


log 


C  ^  P  ^  (C  ^  P)" 

C  -  P  +  m'ol^  >  (C  -  P)" 


[Z(2,  1  -  c  -  p)  -  z(2,  1  -  c  +  p: 

4 


+  Z(2,  1  +  C  -  P)  -  Z(2,  1  +  C  +  P) 


(6.3) 


+  0(a")  } 


where  the  Gamma  Function,  r(x),  is  defined  by 

r  (>r)  =  f  dt 

J  0 

q  is  a  flux  density  given  by 

g=  - 

^  4nkb 

For  the  remaining  larger  values  of  p 


<J)  =  4g  [  — —  Kg(2nna)  cos  (2nnC)  sin(2/2mP) 

^  n  =  l  ^ 

(6.4) 

+  Plog^] 

a 

where  Kq(x)  is  the  modified  Bessel  Function  of  the  second 
kind.  For  relatively  large  values  of  x  (Press,  Flannery, 
Teukolsky,  and  Vetterling  1989), 

K^ix)  ~  -----  e  ^ 

y/2nx 

Thus,  an  implementation  of  Equation  6.4  requires  relatively 


few  terms. 
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Well  in  an  Aquifer 

The  equation  for  the  potential  of  a  well  in  the  simpli¬ 
fied  aquifer  shown  in  Figure  48  can  be  determined  from  the 
solution  of  a  single  well  by  the  method  of  images  (Told 
1959).  Figure  50  illustrates  how  it  works.  The  original  +Q 


Figure  50.  Method  of  Images 


102 


well  can  have  zero  potential  along  the  river  boundary  by 
adding  a  -Q  image  well  reflected  across  the  boundary.  In 
like  manner,  a  +Q  well  can  have  the  impervious  boundary 
condition  preserved  by  adding  a  +Q  image  well  across  the 
boundary.  After  the  original  four  image  wells  have  been 
added,  however,  additional  ones  must  be  added  to  balance  the 
image  wells.  In  fact,  the  exact  solution  is  an  infinite 
number  of  image  wells.  However,  because  the  radius  of 
influence  of  a  well  is  from  500  to  1,000  ft,  a  first  or 
second  order  approximation  is  typically  all  that  is 
required.  However,  the  small  problem  currently  being 
addressed  requires  29  image  wells  for  each  original  well. 
This  gives  15  wells  on  either  side  of  the  river  boundary  to 
provide  zero  potential  at  the  river. 

Equations  6.3  and  6.4  can  now  be  appropriately  applied 
to  the  original  well  and  the  29  image  wells  if  additional 
care  is  taken  in  defining  variables  in  Equation  6.1.  If  a 
well  is  located  at  (Xq,  y^)  and  the  top  of  the  aquifer  is  at 
z^,  then  use 


a  = 

c  = 


1 

2h 

1 

2h 


yjix  - 
-  z) 


iy  - 


Ihe  positions  and  signs  of  their  respective  Q's  for  the 
29  image  wells  are  now  given.  Here,  the  length  of  the  aqui¬ 
fer  is  designated  by  1,  and  the  origin  of  the  coordinate 
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system  as  seen  in  the  plan  view  is  the  lower,  left-hand 
corner  at  the  bottom  of  the  aquifer. 
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Multiple  Wells  in  an  Aquifer 

Multiple  wells  in  an  aquifer  are  simply  handled  by 
summing  the  results  of  the  individual  wells.  Computation¬ 
ally,  it  is  easy  to  specify  a  zero  potential  on  the  part  of 
the  grid  modeling  the  river  by  selecting  a  suitable  value 
for  the  datum.  In  our  simplified  problem,  for  instance,  hp 
=  120  ft.  The  heads  will  all  be  zero  or  negative  and  repre¬ 
sent  the  drawdown  as  a  result  of  the  pumping  of  the  wells. 

Quality  of  Grids  Used 

Usually,  one  can  only  qualitatively  analyze  the  grids 
that  are  used  by  examining  such  things  as  orthogonality, 
aspect  ratio,  smoothness,  etc.,  but,  in  this  case,  since  the 
exact  solution  is  known,  a  quantitative  comparison  can  also 
be  made.  This  offers  an  excellent  opportunity  to  test  some 
of  the  different  options  in  EAGLE  and  different  mappings  to 
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see  what  works  best  for  this  problem.  One  simply  computes 
the  percentage  error  for  each  grid  point  for  each  grid  stud¬ 
ied.  This  gives  an  excellent  way  of  improving  the  grids 
used  in  modeling  wells  in  real-world  problems  rather  than 
simply  using  what  "looks  good." 

Solution  for  a  Partially  Penetrating  Well 

First,  consider  the  solution  for  a  partially  penetrat¬ 
ing  well  in  an  infinite  aquifer  (beyond  the  well's  radius  of 
influence) .  This  will  give  an  indication  as  to  how  to  do 
more  complex  .  oblems  using  wells  as  well  as  show  the  scien¬ 
tific  visa  xization  capabilities  of  FAST  for  this  applica¬ 
tion.  The  problem  solved  is  the  second  well  of  the  simpli¬ 
fy  ed  problem.  The  pertinent  data  are 

Q  =  200  cfm 
k  =  0.1  ft/min 
b  =  80  ft 
=  120  ft 

Well  Radius  =  2  ft 

Radius  of  influence  =  480  ft 

Grid  technique 

O  grid  with  plug.  A  grid  is  constructed  first  by  con¬ 
structing  an  O  type  grid  (21  by  21  by  25)  as  if  doing  a 
fully  penetrating  well  (the  entire  120  ft)  for  one  block  and 
then  for  a  second  block  that  constitutes  a  plug  (40  ft  here 
with  a  6  by  6  by  9  grid)  to  account  for  the  well  being  only 
partially  penetrating.  Figure  51  shows  the  O  grid  repre¬ 
senting  the  soil,  and  Figure  52  shows  the  grid  for  the  plug. 
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Figure  51.  O  Type  Grid 


Note  that  an  O  type  grid  could  have  been  used  for  the  plug 
but  a  different  type  was  chosen. 

Spacing.  The  grids  in  Figures  51  and  52  use  equal 
spacing.  However,  it  is  important  to  use  concentrated 
spacing  toward  the  well  and  where  the  well  bottom  touches 
the  soil.  Tnerefore,  a  second  grid  and  solution  was 
obtained  using  concentrated  spacing.  Figure  53  shows  the 
radial  distribution  of  nodes,  and  Figure  54  shows  the 
distribution  with  depth  of  this  grid. 


Figure  52.  Plug 


Analysis  of  results 

Visualization.  First,  the  results  were  viewed  using 
FAST  to  consider  the  quality  of  the  solution.  It  turns  out 
that  FAST  works  very  well  for  this  purpose.  Figure  55  shows 
a  color  contour  plot  of  total  head  for  a  vertical  section 
(constant  I)  using  the  module  Surfer.  This  plot  is  very 
similar  to  plotting  equipotentials  as  half  of  a  flow  net  in 
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Figure  53.  Radial  Spacing 


2-D  applications.  The  results  show  very  clearly  horizontal 
flow  until  the  water  approaches  the  well.  Figure  56  shows 
some  constant  J  color  contour  plots.  The  constant  color 
bands  indicate  horizontal  flow  as  well.  As  one  approaches 
the  well,  the  constant  J  surfaces  have  different  colors 
showing  other  than  horizontal  flow.  Figure  57  shows  a  constant 


Figure  54.  Depth  Spacing 


Ill 


Figure  55.  Constant  I  Color  Contour  Plot 


Figure  56.  Constant  J  Color  Contour  Plot 


Figure  57.  Constant  K  Color  Contour  Plot 


K  color  contour  plot.  First,  the  problem  is  axisymmetric , 
so  the  solution  must  be  as  well,  which  is  the  case.  Also, 
the  variation  is  approximately  logarithmic  as  should  be. 

Another  type  of  visualization  technique  using  module 
Isolevel  is  shown  in  Figure  58.  The  2-D  problem  has  equipo- 
tential  lines,  but  the  3-D  problem  has  isolevel  surfaces. 
Again,  since  flow  occurs  perpendicular  to  these  surfaces, 
they  will  be  vertical  rings  until  the  well  is  approached. 

Figure  59  shows  the  final  type  of  visualization  used 
(generated  from  module  Tracer) .  It  is  a  particle  trace  for 
a  vertical  section  of  the  soil  near  the  well.  What  the 
aerospace  engineers  use  for  showing  the  flow  of  air  around 
an  aircraft  is  also  an  excellent  tool  for  showing  flow  in 
porous  media.  The  part  of  the  grid  that  represents  the  well 
is  shown  in  red,  and  the  flow  lines  are  shown  in  black. 

Again  horizonal  flow  exists  except  near  the  well. 

The  first  value  of  these  visualization  tools  is  to  aid 
in  deciding  whether  the  solution  appears  reasonable.  If 
incorrect  boundary  conditions  or  poor  answers  result,  this 
becomes  readily  apparent.  A  second  value  is  once  confidence 
has  been  obtained  in  the  model,  specific  details  as  to  what 
is  happening  can  then  be  analyzed. 

Error  analysis.  Equations  6.3  and  6.4  were  used  to 
compute  head  or  potential  for  each  of  the  10,644  nodes,  and 
the  results  were  compared  with  the  computed  results.  The 
break  point  for  a  between  using  Equation  6.3  or  6.4  is 
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Figure  59.  Particle  Trace  Plot 


a  =  0.2 


The  grid  with  equal  spacing  had  a  maximum  percentage  error 
of  24.1  percent.  Since  the  heads  near  r  =  480  ft  approach 
zero,  the  percent  error  was  computed  by 


%  c 


X  100% 


where 

h^  =  computed  head 

hj  =  theoretical  head 

h^  =  maximum  head  at  the  well 

Figure  60  shows  a  color  contour  plot  of  percent  error  for 
this  grid  with  white  being  the  worst  case.  It  is  remarkable 
that  with  even  spacing  the  error  is  very  small  until  you  are 
very  near  the  well.  The  default  spacing  in  EAGLE  was  used, 
and  the  maximum  percent  error  for  this  grid  was  4.6  percent. 


Solution  of  Simplified  Problem 

Now  the  solution  for  the  problem  given  in  Figure  48 
will  be  given.  Because  quality  of  grids  is  an  important 
issue,  different  single  and  multiple  block  configurations 
were  considered.  Because  of  its  potential  to  real-world 
problems  and  the  developed  capability  to  merge  separate 
blocks,  the  selected  solution  is  to  enclose  each  well  in  an 
O  type  box  grid  (Figure  61)  and  merge  them  with  their 
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Figure  60.  Percent  Error  Plot 


Figure  61.  Two-Well  Problem 


respective  plugs  and  each  other.  Two  versions  will  be  pre¬ 
sented  here.  A  large  algebraic  grid  of  21,266  nodes  will  be 
first  considered  and  then  compared  with  a  much  smaller 
elliptic  grid  of  13,266  nodes  using  the  contyp= ' initial ' 
option  in  EAGLE.  Appendix  D  contains  a  description  of  the 
data  required  to  generate  the  elliptic  FEM  grid. 

Analysis  of  results 

Visualization.  First,  the  large  algebraic  grid  at  the 
well  is  shown  in  Figure  62.  The  distribution  of  nodes  from 
the  circular  well  to  the  rectangular  boundary  is  quite  good 
before  smoothing  is  done.  Figures  63  through  65  show  color 
contours  of  total  head  for  J,  I,  and  K  levels,  respectively. 
Here,  the  highest  value  (white)  is  at  the  river,  and  the 
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Figure  62.  Large  Algebraic  Grid  at  Well 
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qure  63.  J  Surface  Color  Contour  Plot 


Figure  64.  I  Surface  Color  Contour  Plot 


Figure  65.  K  Surface  Color  Contour  Plot 


lowest  value  of  head  (black)  is  at  the  wells.  Figure  66 
shows  the  isolevel  plot  for  the  two  well  system.  Figure  67 
shows  flow  lines  (particle  traces)  going  into  the  first 
well.  All  flow  starts  at  the  river  and  ends  at  one  of  the 
two  wells.  One  surprising  result  is  that  the  left  most  flow 
line  skips  the  smaller  well  and  goes  for  the  bigger  one. 

Error  analysis.  The  large  algebraic  grid  of  21,266 
nodes  created  from  essentially  putting  together  two  single 
well  grids  had  a  maximum  percentage  error  of  5.3  percent. 
Actually,  some  of  this  error  is  attributed  to  the  truncation 
in  the  number  of  image  wells  and  some  is  due  to  numerical 
imperfections.  It  seemed  plausible  that  this  grid  could  be 
significantly  reduced  in  size  by  applying  the  elliptic  grid 
generation  techniques  of  EAGLE  to  get  close  to  the  same 
result.  Figure  68  shows  a  K  surface  for  the  first  well  for 
the  smaller  elliptic  grid  of  13,266  nodes.  This  grid 
yielded  a  maximum  percentage  error  of  5.7  percent  which 
substantiates  the  hypothesis. 

Real-World  Example 

A  real-world  example  showing  the  solution  process  will 
now  be  presented. 

Description  of  Problem 

The  problem  shown  in  plan  view  in  Figure  69  consists  of 
a  part  of  an  aquifer  containing  a  river  crossing  through  the 
region  with  partially  penetrating  wells  pumping  water  out  of 
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Figure  66.  Isolevel  Plot 


Figure  67.  Flow  Lines  Plot  for  One  Well 


Figure  68.  Elliptic  Grid 

the  system.  Line  segment  EF  is  an  impervious  wall  such  as  a 
slurry  trench  found  at  certain  sensitive  military  arsenals 
and  construction  sites.  The  region  with  the  two  wells  is 
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Figure  69.  Plan  View  of  Aquifer 
highly  anisotropic  and  has  the  following  permeability  data 
(see  Figure  70) ; 

0e  =  21° 

=  -53° 

=  10° 

=0.01  /c/min 
kp2  =  0  •  015  f  t/min 


k^,  =  0.001  f  t/min 


z 


Figure  70.  Permeability  Orientation 
The  region  under  the  river  is  a  rather  impervious  clay  hav 
ing  the  soil  properties: 

0,  =  0° 

=  0° 

ilr^  = 

=  0.00001  ft/min 
kp2  =  0.00001  ft/min 
kpj  =  0.00001  /t/min 

Finally,  the  region  with  the  slurry  trench  has  soil  proper 
ties  of  a  pervious  sand: 
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e.  =  0“ 

4‘e  =  0“ 

'ife  =  0“ 

kpj^  =  0.1  f  t/min 
kp2  =  0.1  f  t/min 
kpj  =  0.1  f  t/min 

The  thickness  of  the  aquifer  varies  between  approximately 
200  to  350  ft. 

FEM  Grid 

This  problem  provides  an  excellent  opportunity  to  test 
the  techniques  and  concepts  developed  thus  far.  The  region 
was  divided  into  16  blocks  or  subregions  as  shown  in 
Figure  71.  Thirteen  subregions  are  visible  with  three  plugs 
making  partially  penetrating  wells  not  shown.  Decoupling 
the  geometry  makes  it  much  easier  to  generate  the  grid  at 
times.  It  is  certainly  not  necessary  to  have  so  many  blocks 
as  the  mapping  capability  of  EAGLE  is  extensive.  However, 
it  was  done  this  way  to  give  the  approach  a  thorough  test. 
Once  the  data  files  for  the  surf  and  grid  modules  of  EAGLE 
were  developed  for  a  particular  type  of  region,  it  was  very 
easy  to  modify  these  files  for  the  next  similar  type  region. 
After  all  the  pieces  were  generated,  it  was  a  simple  matter 
to  combine  the  pieces  into  a  single,  complete  FEM  grid  using 
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Figure  71.  Subregions 


the  program  developed  as  part  of  this  work.  Figure  72  shows 
a  plan  view  of  the  grid,  and  Figure  73  shows  a  perspective 
view  of  the  grid  and  aquifer  system.  Eleven  layers  (K 
levels)  were  used  with  the  resulting  FE  grid  having 
11,578  nodes  and  9,855  elements. 

With  the  head  at  the  river  being  300  ft,  the  wells  each 
having  a  penetration  of  100  ft,  and  a  head  of  260  ft  at  the 
well,  the  following  data  was  used  to  generate  the  FE  grid 
for  the  problem  assuming  the  16  pieces  have  already  been 
generated : 


com 

16 

aquif 1 . egl 
1  21.  -53.  10. 
aquif2 .egl 
1  21.  -53.  10. 
aquif3 .egl 
1  21.  -53.  10. 
aquif 4 . egl 
1  21.  -53.  10. 
aquif 5. egl 

1  21.  -53.  10. 
aquife.egl 

2  0  0  0 
aquif? . egl 
3  0  0  0 
aquif 8 . egl 
3  0  0  0 
aquifO .egl 
3  0  0  0 
aquif 10 . egl 
3  0  0  0 
aquif 11 . egl 
3  0  0  0 
aquif 12 . egl 
3  0  0  0 
aquif 13 .egl 
3  0  0  0 
aquif 14 .egl 

1  21.  -53.  10. 
aquif 15 . egl 

1  21.  -53.  10. 
aquif 16 . egl 

3  0  0  0 

ban 

be 

2116  23  181  260. 
be 

2  1  1  9  23  1  9  12  260. 
be 

2  1  1  10  23  1  11  2  0. 
be 

5116  23  181  260. 
be 

5  1  1  9  23  1  9  12  260. 
be 

5  1  1  10  23  1  11  2  0. 
be 

8116  21  181  260. 
be 

8  1  1  9  21  1  9  12  260. 
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be 

8  1  1  10  21  1  11  2  0. 
be 

6  1  1  11  23  5  11  1  300 
out 

aquifer 

Aquifer  with  Uneonfined  Flow 

0. 

.01  .015  .001 

.0001  .0001  .0001 

.1  .1  .1 

end 


The  grid  around  eaeh  well  was  generated  by  going 
direetly  from  a  eirele  to  a  reetangular  type  region  with 
eare  to  use  gradually  inereasing  spaeing  (Figure  74) .  For 
some  applieations  there  may  be  too  mueh  skewness,  and  an 
alternate  plan  sueh  as  a  five-region  system  must  be  imple¬ 
mented.  However,  as  will  be  shown  later,  good  results  were 
achieved  in  this  application  using  the  original  approach. 

The  impervious  wall  can  be  easily  handled  as  shown  in 
Figure  75.  All  that  must  be  done  is  to  have  nodes  at  dif¬ 
ferent  positions  on  the  sides  of  the  wall. 


Presentation  of  Results 

The  problem  was  first  run  as  having  confined  flow  and  a 
homogeneous,  isotropic  medium.  Figure  76  shows  a  color 
contour  plot  of  potential  for  the  top  value  of  k  =  11  with 
white  being  full  potential  and  black  representing  the  lowest 
value  of  potential.  An  exception  is  the  river  which  is  at 
full  potential  but  is  painted  blue  for  visualization  effec¬ 
tiveness.  Radial  flow  at  the  well  is  seen  which  is  to  be 
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Figure  75.  Grid  at  Impervious  Wall 


expected.  It  is  also  interesting  to  see  the  effect  of  the 

impervious  wall  on  the  head  distribution. 

The  unconfined  flow  problem  was  then  run  for  the  real 

soil  conditions,  requiring  a  running  time  of  3763.2  sec  on 

the  Cray  YMP  and  five  iterations  for  convergence. 
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Figure  76.  Homogeneous,  Isotropic  Medium 


An  important  aspect  of  groundwater  flow  is  the  amount  of 
flow  from  the  wells.  The  three  wells  pump  a  total  of 
512.7  cfm,  which  is  within  their  range  of  capacity. 

Figure  77  shows  the  color  contour  plot  for  unconfined  flow. 
Now,  the  less  pervious  region  with  the  two  wells  has  more 
head  loss  which  is  correct.  The  free  surface  algorithms  can 
be  tested  by  viewing  the  free  surface  at  a  well.  Figure  78 
shows  a  portion  of  the  grid  near  one  of  the  wells.  The 
results,  including  the  exit  line,  look  as  they  should.  Many 
other  plots  can  be  obtained,  but  these  are  left  to  the 
interested  user  of  the  program  to  generate. 
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Figure  77.  Actual  Medium 


CHAPTER  VII 


SUMMARY  AND  CONCLUSION 

The  techniques  and  tools  ordinarily  used  by  aerospace 
engineers  were  successfully  applied  to  2-D  and  3-D  seepage 
and  groundwater  modeling. 

2-D  Flow  Nets 

First,  the  techniques  used  to  generate  an  orthogonal 
grid  were  broadened  in  an  innovative  way  to  produce  computer 
generated  flow  nets  of  superior  quality  to  those  from  other 
techniques.  When  compared  with  less  theoretically  based 
algorithms,  the  results  using  the  author's  work  proved  to  be 
significantly  more  accurate. 

3-D  Modeling 

Next,  computational  techniques  using  the  FEM  were 
developed  for  3-D  flow  applications  to  apply  correct  bound¬ 
ary  conditions  for  unconfined  flow  and  dodge  the  problems 
that  occur  at  exit  lines.  The  avoidance  of  a  complicated 
2-D  surface  extrapolation  for  the  exit  line  proved 
especially  helpful.  This  allowed  the  successful  completion 
of  a  3-D  seepage  and  groundwater  model  based  on  the  interna¬ 
tionally  known  2-D  seepage  package  developed  by  this  author. 
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A  very  innovative  technique  and  program  were  next 
developed  to  allow  the  user  to  first  generate  pieces  of 
structured  grid  using  EAGLE  and  then  combine  the  pieces  in 
an  efficient  manner  to  produce  the  unstructured  FEM  grid. 
This  has  several  advantages: 

1.  It  makes  it  easy  to  apply  boundary  conditions. 

2.  Decoupling  the  geometry  makes  it  easier  and  faster 
to  generate  the  resulting  FEM  grid. 

3.  Boundaries  between  blocks  only  have  to  match.  Cuts 
and  differences  in  computational  coordinates  are 
automatically  fixed. 

4.  The  grids  are  very  aesthetic,  which  is  important 
when  showing  work  to  upper  management. 

Finally,  the  structure  of  the  original  pieces  is  pre¬ 
served  and  output  is  placed  in  FAST  format  for  later  visual 
ization.  The  innovative  aspect  of  this  work  is  that  for 
unconfined  flow  problems  the  grid  collapses  in  an  arbitrary 
manner  to  the  free  surface,  and  the  structure  is  difficult 
to  maintain.  In  fact,  the  only  thing  that  is  used  is  the 
number  of  blocks  and  the  dimensions  of  each  block  to  recon¬ 
struct  both  the  collapsed  grid  and  results  files. 

All  the  above  tools  were  then  applied  to  the  16-block 
problem  of  groundwater  flow  in  an  aquifer.  It  was  gratify¬ 
ing  to  see  that  the  entire  process  really  works! 
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This  appendix  contains  the  listing  and  description  of 
three  subroutines  which  are  a  FORTRAN  implementation  of  the 
free  surface  algorithms  developed  as  part  of  this  work. 


Subroutine  FREES 

Subroutine  FREES  determines  the  final  position  of  the 
free  surface,  and  its  FORTRAN  listing  follows. 


SUBROUTINE  FREES 
C 
C 

C  THIS  SUBROUTINE  COMPUTES  THE  POSITION  OF  THE  FREE 

C  SURFACE . 

C 

C 

C  DESCRIPTION  OF  NBC  VALUES  OF  100  OR  GREATER. 

C 

C  100  -  ABOVE  THE  FS  WITH  FLOW  NOT  GIVEN 

C  200  -  ABOVE  THE  FS  WITH  FLOW  GIVEN 

C  300  -  ON  FS  WITH  FLOW  NOT  GIVEN 

C  400  -  ON  FS  WITH  FLOW  GIVEN 

C  500  -  BELOW  FS  WITH  FLOW  NOT  GIVEN 

C  600  -  BELOW  FS  WITH  FLOW  GIVEN 

C 
C 

COMMON  /  GRIDl  /  NUMNP,  NUMEL,  NUMMAT,  X(IOOO), 

&  Y(IOOO),  Z(IOOO),  FX(IOOO),  NBC(IOOO),  FLOW(IOOO), 

&  HLAST(IOOO) 

COMMON  /  GRID2  /  XK1(12),  XK2(12),  XK3(12), 

&  NP(9,  900),  THETA(900),  PHI{900),  PSI(900) 

COMMON  /  FREE  /  FRX(IOOO),  FRY (1000),  FRZ(IOOO), 

&  COUNT (1000) 

COMMON  /  BANARG  /  MBAND,  NUMBLK,  R(IOOO),  C(120,60), 
&  ND,  ND2 

DIMENSION  IADJ(8,  3) 

DATA  lADJ  /2,  3,  4,  1,  6,  7,  8,  5,  4,  1,  2,  3,  8,  5, 

&  6,  7,  5,  6,  7,  8,  1,  2,  3,  4/ 

INITIALIZE  DATA. 

DO  100  1=1,  NUMNP 
COUNT (I)  =  1.E30 
FRX(I)  =  0. 

FRY (I)  =  0. 

FRZ(I)  =  0. 

100  CONTINUE 
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DETERMINE  THE  FREE  SURFACE  NODES. 

DO  190  N  =  1,  NUMEL 

DO  180  I  =  1,  8 

N1  =  NP(I,  N) 

II  =  IADJ(I,  1) 

N2  =  NP(II,  N) 

II  =  IADJ(I,  2) 

N3  =  NP(II,  N) 

II  =  IADJ(I,  3) 

N4  =  NP(II,  N) 

PHI  =  R(N1)  -  Z(N1) 

IF  (PHI)  120,  110,  180 

FIX  ALL  BUT  NBC  =  2  NODE. 

110  IF  (NBC(Nl) .EQ.2)  GO  TO  180 
FRX(Nl)  =  X(N1) 

FRY(Nl)  =  Y(N1) 

FRZ(Nl)  =  Z(N1) 

COUNT (Nl)  =  0. 

GO  TO  180 

CONSIDER  LINE  SEGMENT  N1-N2. 

120  CALL  FSPT(N1,  N2) 

CONSIDER  LINE  SEGMENT  N1-N3 . 

CALL  FSPT(N1,  N3) 

CONSIDER  LINE  SEGMENT  N1-N4 . 

CALL  FSPT(N1,  N4) 

180  CONTINUE 
190  CONTINUE 

SET  FLAGS  AND  MOVE  FREE  SURFACE  NODES. 
DO  300  1=1,  NUMNP 
NBI  =  NBC(I) 

IF  (COUNT(I) .LT. 1.E29)  THEN 
NBI  =  300 

IF  (NBC(I) .NE.O)  NBI  =  400 
COUNT{I)  =  1. 
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X(I)  =  FRX(I) 

Y(I)  =  FRY (I) 

Z(I)  =  FRZ(I) 

R(I)  =  Z(I) 

C 

C  FIX  FLOW  FOR  NBC  =  -1. 

C 

IF  (NBC(I) .EQ.-l)  FLOW(I)  =  FX(I) 

ELSE 

COUNT ( I )  =  0 . 

IF  (R(I) .LT.Z(I) )  THEN 
NBI  =  100 

IF  (NBC(I) .NE.O)  NBI  =  200 
ELSE 

NBI  =  500 

IF  (NBC(I) .NE.O)  NBI  =  600 
ENDIF 
ENDIF 
C 

NBC (I)  =  NBI 
C 

300  CONTINUE 
C 

C  FIX  ELEMENTS  CLOSE  TO  THE  FREE  SURFACE. 

C 

CALL  FIXEL 
C 

RETURN 

END 

The  codes  100-600  in  the  comments  are  flags  used  to 
determine  how  the  node  is  to  be  considered  upon  print-out  of 
results.  A  node  is  below,  on,  or  above  the  free  surface 
with  flow  given  or  not  given.  Variables  FRX,  FRY,  and  FRZ 
contain  the  coordinates  of  the  new  position  of  a  free  sur¬ 
face  node,  and  the  value  of  COUNT  determines  if  the  node  is 
a  free  surface  node.  COUNT  is  initialized  to  10^°,  and  if 
it  is  set  to  a  smaller  value,  the  node  is  a  free  surface 
node.  Do  loops  180  and  190  check  each  line  segment  of  each 
element  to  find  a  line  segment  where  the  pressure  head  is 
less  than  zero  on  the  first  node  and  greater  than  zero  on 
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the  second  node.  When  such  a  line  segment  is  found,  subrou¬ 
tine  FSPT  is  called  to  compute  the  new  position  of  the  free 
surface  node. 


The  next  section  (do  loop  300)  does  three  things: 

1.  Changes  the  value  of  COUNT  to  be  1  for  a  free 
surface  node  and  0  otherwise. 

2.  Puts  in  the  array  NBC  the  print  codes  (100-600) 
describing  the  status  of  each  node. 

3-  Changes  the  (x,  y,  z)  coordinates  of  the  free 

surface  nodes  to  their  respective  new  free  surface 
values . 

Finally,  subroutine  FIXEL  is  called  to  fix  the  elements 
to  be  either  all  below  or  all  above  the  free  surface. 


Subroutine  FSPT 

Subroutine  FSPT  looks  for  a  free  surface  point  between 
nodes  N1  and  N2.  Its  listing  is  now  given. 


SUBROUTINE  FSPT(N1,  N2) 

C 

C 

C  THIS  SUBROUTINE  CHECKS  LINE  SEGMENT  M1-N2  FOR  A 

C  POTENTIAL  FREE  SURFACE  NODE.  IF  SO,  THE  NEW 

C  POSITION  OF  NODE  N1  IS  COMPUTED. 

C 

C 

COMMON  /  GRIDl  /  NUMNP,  NUMEL,  NUMMAT,  X(IOOO), 

&  Y(IOOO),  Z(IOOO),  FX(IOOO),  NBC(IOOO),  FLOW(IOOO), 

&  HLAST(IOOO) 

COMMON  /  GRID2  /  XK1(12),  XK2(12),  XK3(12), 

&  NP(9,  900),  THETA(900),  PHI (900) ,  PSI(900) 

COMMON  /  FREE  /  FRX(IOOO),  FRY(IOOO),  FRZ(IOOO), 

&  COUNT(IOOO) 

COMMON  /  BANARG  /  MBAND,  NUMBLK,  R(IOOO),  C(120,60), 
&  ND,  ND2 
C 

C  SET  EXIT  POINT  PARAMETER. 

C 

BETA  =  . 8 

C 

i72 
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PHI  =  R(N1)  -  Z(N1) 

PH2  =  R(N2)  -  Z(N2) 

C 

IF  (PH2)  140,  100,  110 

COMPUTE  S  USING  THE  SPECIAL  EXIT  LINE  COMPUTATION. 

100  IF  ( (NBC{N1)  .NE.2)  .OR.  (NBC(N2) -NE.2) )  GO  TO  140 
DZZ  =  Z(N2)  -  Z(N1) 

IF  (DZZ.EQ.O.)  THEN 
S  =  1. 

ELSE 

SO  =  PHI  /  DZZ 
S  =  (1  -  BETA)  *  SO  +  BETA 
ENDIF 
GO  TO  120 

DON'T  LET  A  SURFACE  OF  SEEPAGE  NODE  BE  MOVED  OFF 
THE  SURFACE  OF  SEEPAGE. 

110  IF  (NBC(Nl) .EQ. 2)  GO  TO  140 

COMPUTE  S  USING  THE  FREE  SURFACE  NODE  COMPUTATION. 

S  =  PHI  /  (PHI  -  PH2) 

120  IF  ( (S.LT. 0. ) .OR. (S.GT. 1. ) )  S  =  1. 

DX  =  (X(N2)  -  X(N1) )  *  S 

DY  =  (Y(N2)  -  Y(N1) )  *  S 

DZ  =  (Z  (N2)  -  Z  (Nl)  )  *  S 

RR  =  DX  *  DX  +  DY  *  DY  +  DZ  *  DZ 

KEEP  THE  ONE  WITH  THE  SMALLEST  MOVEMENT. 

IF  (RR.GE.COUNT(Nl) )  GO  TO  140 
FRX(Nl)  =  X(N1)  +  DX 

FRY(Nl)  =  Y(N1)  +  DY 

FRZ(Nl)  =  Z(N1)  +  DZ 

COUNT (Nl)  =  RR 
C 

140  RETURN 
END 


The  basic  idea  of  this  subroutine  is  to  compute  the 
value  of  the  parameter  s,  where 


0  s  s  1 
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such  that  the  pressure  head  is  zero.  That  is,  use  Equation 
4.28  for  a  node  not  on  the  exit  face  or  Equations  4.29  and 
4.30  for  a  node  on  the  exit  face  (NBC  =  2).  The  chanqes  in 
coordinate  values  (variables  DX,  DY,  and  DZ)  required  to 
move  point  N1  to  the  new  free  surface  position  are  then 
computed.  Distance  squared  (in  variable  RR)  is  then 
computed  from  these  values  and  compared  with  the  current 
value  of  COUNT.  If  RR  is  less  than  COUNT,  the  position  of 
the  free  surface  is  computed  along  this  line  segment  and 
stored  in  FRX,  FRY,  and  FRZ.  COUNT  is  then  updated  to  RR. 
This  position  of  the  free  surface  is  kept  for  node  N1  until 
a  smaller  RR  comes  along  as  a  result  of  considering  other 
line  segments. 


Subroutine  FIXEL 

This  subroutine  performs  the  iterative  process  of  rede¬ 
fining  the  coordinates  of  the  node  points  so  elements  are 
either  all  below  the  free  surface  or  collapsed  onto  it.  Its 
listing  is  given. 


SUBROUTINE  FIXEL 
C 
C 

C  THIS  SUBROUTINE  FIXES  THE  ELEMENTS  CROSSED  BY  THE 

C  FREE  SURFACE  SO  AN  ELEMENT  IS  EITHER  ALL  BELOW  THE 

C  FREE  SURFACE  OR  COLLAPSED  ONTO  IT. 

C 

C 

COMMON  /  GPIDl  /  NUMNP,  NUMEL,  NUMMAT,  X(IOOO), 

&  Y(IOOO),  Z(IOOO),  FX(IOOO),  NBC(IOOO),  FLOW(IOOO), 

&  HLAST(IOOO) 

COMMON  /  GRID2  /  XK1(12),  XK2(12),  XK3(12),  NP(9, 

&  900),  THETA(900),  PHI(900),  PSI(900) 
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COMMON  /  FREE  /  FRX(IOOO),  FRY (1000),  FRZ(IOOO), 

&  COUNT (1000) 

COMMON  /  BANARG  /  MBAND,  NUMBLK,  R(IOOO),  C(120,60), 
&  ND,  ND2 

DIMENSION  ISEG(2,  12) 

DATA  ISEG  /I,  2,  2,  3,  3,  4,  4,  1,  1,  5,  2,  6,  3,  7, 

&  4,  8,  5,  6,  6,  7,  7,  8,  8,  5/ 

KEEP  ITERATING  UNTIL  THERE  IS  NO  CHANGE. 

ICOUNT  =  0 


100  ICHANG  =  0 

ICOUNT  =  ICOUNT  +  1 
IF  ( ICOUNT. GT. 200)  THEN 

PRINT*,  '  AFTER  200  ITERATIONS  THE  SUBROUTINE', 

&  '  FIXED  ALGORITHM  HAD  NOT  CONVERGED. ' 

RETURN 
ENDIF 

DO  600  N  =  1,  NUMEL 

CHECK  IF  THE  FREE  SURFACE  INTERSECTS  THE  ELEMENT. 


PMIN  =  1.E30 
PMAX  =  -  PMIN 
DO  200  I  =  1,  8 
II  =  NP(I,  N) 

P  =  R(II)  -  Z(II) 

PMIN  =  AMIN1(PMIN,  P) 
PMAX  =  AMAX1(PMAX,  P) 
200  CONTINUE 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


IF  THE  ELEMENT  IS  ABOVE  THE  FREE  SURFACE  AND  DOES 
NOT  TOUCH  IT,  DON'T  BOTHER  WITH  IT. 

IF  (PMAX.LT.O.)  GO  TO  600 

IF  THE  ELEMENT  IS  BELOW  THE  FREE  SURFACE,  DON'T 
BOTHER  WITH  IT. 

IF  (PMIN.GE.O.)  GO  TO  600 

THE  ELEMENT  IS  CROSSED  BY  THE  FREE  SURFACE,  SO 
CHECK  EACH  LINE  SEGMENT. 

DO  500  J  =  1,  12 

II  =  ISEG(1,  J) 

N1  =  NP(II,  N) 

II  =  ISEG(2,  J) 
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=  NP(II,  N) 

IF  LINE  SEGMENT  N1-N2  IS  A  FREE  SURFACE  NODE  AND  A 
NODE  ABOVE  THE  FREE  SURFACE,  COLLAPSE  THE  NODE 
ABOVE  THE  FREE  SURFACE  TO  THE  EXIT  POINT  NODE. 

300  K  =  1,  2 

CONSIDER  BOTH  N1-N2  AND  N2-N1, 

(K.EQ.2)  THEN 
NN  =  N1 
N1  =  N2 
N2  =  NN 
ENDIF 

IF  ( (COUNT(Nl) .GT.O.) .AND. (R(N2) .LT.Z(N2) ) )  THEN 
X(N2)  =  X(N1) 

Y(N2)  =  Y(N1) 

Z(N2)  =  Z(N1) 

R(N2)  =  R(N1) 

COUNT (N2)  =  2. 

FLAG  THAT  A  CHANGE  HAS  BEEN  MADE. 

I CHANG  =  1 
ENDIF 

300  CONTINUE 
500  CONTINUE 
600  CONTINUE 

IF  A  CHANGE  WAS  MADE,  REPEAT  THE  PROCESS. 

IF  (ICHANG.EQ. 1)  GO  TO  100 

RETURN 

END 


C 

C 

C 

C 

C 


N2 


DO 


IF 


First,  note  that  a  maximum  of  200  iterations  are 
allowed.  Next,  the  maximum  and  minimum  pressure  heads  for 
the  element  are  computed.  If  they  are  all  greater  than  or 
equal  to  zero,  the  element  is  completely  under  the  free 
surface.  If  they  are  all  less  than  zero,  the  element  is 
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completely  above  the  free  surface  and  does  not  touch  it.  In 
each  of  these  cases  the  element  needs  no  further  action. 
Otherwise,  the  algorithm  given  at  the  end  of  Chapter  IV  is 
executed. 

Here,  each  line  segment  of  the  element  is  considered 
for  the  first  node  being  a  free  surface  node  and  the  second 
node  of  the  line  segment  being  above  the  free  surface 
(pressure  head  less  than  zero) .  When  such  a  line  segment  is 
found,  node  2  is  collapsed  to  node  1  of  the  line  segment, 
and  node  2  is  flagged  as  a  special  free  surface  node  (COUNT 
=  2)  . 

After  going  through  all  the  elements  (do  loop  600)  with 
no  further  collapses  occurring,  the  iteration  process  is 
terminated . 
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This  appendix  documents  the  program  written  to  combine 
blocks  generated  by  the  EAGLE  program  into  a  consolidated 
FEM  grid.  This  approach  makes  it  very  easy  to  generate 
separate  pieces  and  then  combine  them  later.  This  approach 
also  makes  it  significantly  easier  than  usual  to  apply 
boundary  conditions. 


MAIN  Program 

The  MAIN  program  is  the  driver  for  the  rest  of  the 
program  and  is  given: 


C  THIS  PROGRAM  COMBINES  EAGLE  BLOCKS  INTO  ONE  BIG 

C  FEM  GRID. 

C 

C 

PARAMETER  (ND  =  1000) 

COMMON  /OUTPUT/  X(ND),  Y(ND),  Z(ND),  IBC(ND),  BV(ND), 

&  IX(9,  ND) ,  THETA(ND),  PHI(ND),  PSI(ND),  NUMNP,  NUMEL, 
&  NUMMAT,  NODVEL 
CHARACTER  COM  *  3,  ANAME  *  20 
C 

NUMNP  =  0 
NUMEL  =  0 
NUMMAT  =  0 
NODVEL  =  0 
C 

OPEN  (7,  ACCESS^ ' DIRECT' ,  RECL=16) 

OPEN  (8,  FORM= 'UNFORMATTED' ) 

OPEN  (9,  FORM= 'UNFORMATTED' ) 

C 

C  PROCESS  THE  COMMANDS 

C 

100  PRINT*,  '  ' 

PRINT*,  'COMMAND?' 

READ  (*,  110)  COM 

110  FORMAT (A) 


IF 

( (COM. 

,EQ. 

'  INP' 

) 

.OR. 

(COM. 

EQ. 

'  inp ' 

)  ) 

GO 

TO 

400 

IF 

( ( COM . 

■  eq. 

'COM' 

) 

.OR. 

(COM. 

.  EQ. 

'  com ' 

)  ) 

GO 

TO 

200 

IF 

( ( COM . 

.  EQ. 

'OUT' 

) 

.OR. 

(COM. 

EQ. 

’  out ' 

)  ' 

GO 

TO 

500 

IF 

( ( COM . 

.EQ. 

'BAN' 

) 

.OR. 

(COM. 

,EQ. 

'  ban ' 

)  ) 

GO 

TO 

300 

IF 

( (COM. 

.EQ. 

'BC  ' 

) 

.OR. 

(COM. 

,  EQ. 

'be  ' 

)  ) 

GO 

TO 

600 

IF 

( (COM. 

.  EQ. 

'CLE' 

) 

.OR. 

(COM. 

,  EQ. 

'  cle ' 

)  ) 

GO 

TO 

700 

IF 

( (COM. 

.  eq. 

'  END' 

) 

.OR. 

(COM. 

•  EQ. 

'  end ' 

)  ) 

GO 

TO 

800 
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PRINT  120 


120  F0RMAT(/  '  INP 

&  '  COM 

&  '  OUT 

&  '  BAN 

&  '  BC 

&  '  CLE 

&  '  END 


GO  TO  100 


INPUT  A  SINGLE  FILE.'  / 
COMBINE  SEVERAL  FILES.'  / 
OUTPUT  A  FILE. '  / 

BANDWIDTH  MINIMIZE.'  / 

APPLY  BOUNDARY  CONDITIONS. ’  / 
CLEAR  PREVIOUS  WORK. '  / 

END  EXECUTION  OF  PROGRAM.') 


COMBINE  THE  GRIDS. 


200  PRINT*,  '  ' 

PRINT*,  'NUMBER  OF  FILES' 
READ  (*,  *)  NOFILE 
CALL  COMBIN (NOFILE,  IBL) 
GO  TO  100 


BANDWIDTH  MINIMIZE  THE  GRID. 


300  CALL  ADJAC 
CALL  BANMIN 
GO  TO  100 


READ  A  GRID. 


400  CALL  COMBIN(l,  IBL) 
GO  TO  100 


OUTPUT  THE  GRID. 

500  PRINT*,  '  ' 

PRINT*,  'OUTPUT  FILE  NAME  FOR  SEEPAGE/GROUNDWATER', 

&  '  MODEL  (WITHOUT  EXTENSIONS)?’ 

READ  (*,  110)  ANAME 

NN  -  NONBLK (ANAME) 

OPEN  (2,  FILE=ANAME(1  :  NN)  //  '.sep', 

&  STATUS^ ' UNKNOWN ' ) 

OPEN  (3,  FILE=ANAME(1  :  NN)  //  '.dim', 

&  STATUS^ ' UNKNOWN ' ) 

CALL  OUTFEM(IBL) 

CLOSE  (2) 

CLOSE  (3) 

GO  TO  100 

APPLY  BOUNDARY  CONDITIONS. 

600  PRINT*,  '  ' 

PRINT*,  'BLOCK  NUMBER,  IIB,  JIB,  KIB,  I2B,  J2B,  K2B,', 
&  '  IBC,  BV?' 

READ  (*,  *)  IBLNO,  IIB,  JIB,  KIB,  I2B,  J2B,  K2B,  NBC, 

&  BCVAL 
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CALL  BC(IBLNO,  IIB,  JIB,  KIB,  I2B,  J2B,  K2B,  NBC, 
&  BCVAL) 

GO  TO  100 

CLEAR  PREVIOUS  WORK. 

700  NUMNP  =  0 
NUMEL  =  0 
NUMMAT  =  0 
NODVEL  =  0 
REWIND  7 
REWIND  8 
REWIND  9 
GO  TO  100 

TERMINATE  PROCESSING. 

800  CALL  DONE 
STOP 
END 


There  are  seven  commands  provided  to  accomplish  the 
task  of  creating  a  single  FEM  grid  from  individual  blocks  as 
follows : 

1.  Input  -  Input  a  single  EAGLE  grid  file. 

2.  Combine  -  Combine  a  group  of  EAGLE  grid  files. 

3.  Output  -  Output  the  FEM  grid. 

4.  Bandwidth  -  Minimize  the  bandwidth  of  the  FEM  grid. 

5.  BC  -  Apply  boundary  conditions, 

6.  Clear  -  Clear  previous  work. 

7.  End  -  End  execution  of  program. 

A  description  of  the  subroutines  supporting  these  commands 
is  now  given. 


Subroutine  ADJAC 

This  subroutine  produces  an  adjacency  table  giving  the 
nodes  that  are  adjacent  to  a  given  node.  It  is  used  in  the 
bandwidth  minimization  process  and  is  as  follows: 
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SUBROUTINE  ADJAC 


C 

C 

C  THIS  SUBROUTINE  COMPUTES  AN  ADJACENCY  TABLE  FOR 

C  THE  GENERATED  MESH. 

C 

C 

PARAMETER  (ND  =  1000) 

COMMON  /OUTPUT/  X(ND),  Y(ND),  Z(ND),  IBC(ND),  BV(ND), 

&  IX(9,  ND) ,  THETA(ND),  PHI(ND),  PSI(ND),  NUMNP,  NUMEL, 
&  NUMMAT,  NODVEL 
COMMON  /ADJ/  JZ(ND,  9) 

DIMENSION  IADJ(8,  3) 

DATA  lADJ  /2,  3,  4,  1,  6,  7,  8,  5,  4,  1,  2,  3,  8,  5, 

&  6,  7,  5,  6,  7,  8,  1,  2,  3,  4/ 

C 

C 

C  JZ  -  ARRAY  CONTAINING  THE  NEIGHBORING  NODES  FOR 

C  EACH  NODE  (ADJACENCY  TABLE) .  NO  MORE  THAN 

C  EIGHT  NEAREST  NEIGHBORS  MAY  EXIST  FOR  A 

C  GIVEN  NODE. 

C 

C 

C  ZERO  THE  JZ  ARRAY. 

C 

DO  200  N  =  1,  NUMNP 

JZ(N,  9)  =  0 

200  CONTINUE 

CONSIDER  EACH  ELEMENT. 

DO  1300  N  =  1,  NUMEL 

CHECK  EACH  NODE  OF  EACH  ELEMENT. 

DO  1200  L  ==  1,  8 
C 

NODP  =  IX(L,  N) 

C 

DO  1100  M  =  1,  3 
C 

LC  =  IADJ(L,  M) 

NODS  -  IX (LC,  N) 

KK  -  JZ(NODP,  9) 

C 

C  PLACE  THE  SECONDARY  NODE  INTO  THE  ROW  OF  JZ 

C  ORRESPONDING  TO  THE  PRIMARY  NODE,  IF  IT  HASN'T 

C  ALREADY  BEEN  DONE. 

C 

IF  (KK)  700,  700,  900 
C 

700  JZ(NODP,  1)  =  NODS 
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JZ(NODP,  9)  =  1 
GO  TO  1100 
C 

900  DO  1000  K  =  1,  KK 

IF  (JZ(NODP,  K).EQ.NODS)  GO  TO  1100 
1000  CONTINUE 

K  =  KK  +  1 
IF  (K.GT.8)  THEN 

PRINT*,  'THE  NODE',  NODP,  '  HAS  MORE  THAN  EIGHT', 
&  '  LINES  GOING  INTO  IT. ' 

CALL  DONE 
STOP 
ENDIF 

1050  JZ(NODP,  K)  =  NODS 
JZ(NO.DP,  9)  =  K 
C 

1100  CONTINUE 
C 

1200  CONTINUE 
C 

1300  CONTINUE 
C 

RETURN 

END 


Subroutine  BANMIN 

This  subroutine  does  the  bandwidth  minimization  for  the 
FEM  grid.  Its  listing  is  now  given: 


SUBROUTINE  BANMIN 
C 
C 

C  THIS  SUBROUTINE  REDUCES  THE  BANDWIDTH  OF  A  FEM 

C  CONFIGURATION. 

C 

C 

PARAMETER  (ND  =  1000) 

COMMON  /OUTPUT/  X(ND),  Y(ND),  Z(ND),  IBC(ND),  BV(ND), 

&  IX(9,  ND) ,  THETA(ND),  PHI(ND),  PSI (ND) ,  NUMNP,  NUMEL, 
&  NUMMAT,  NODVEL 
COMMON  /ADJ/  JZ(ND,  9) 

COMMON  /SCRAT/  NODE(ND),  NNS(ND),  NOS(ND) 

C 

C 

C  NOS  -  ARRAY  CONTAINING  THE  NEW  NODE  NUMBERS  WITH 

C  THE  OLD  NODE  NUMBERS  AS  SUBSCRIPTS. 

C 
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NNS 


ARRAY  CONTAINING  THE  OLD  NODE  NUMBERS  WITH 
THE  NEW  NODE  NUMBERS  AS  SUBSCRIPTS. 


NODE  -  ARRAY  CONTAINING  THE  BEST  NOS  ARRAY. 
NBWMIN  -  MINIMUM  MAXIMUM  DIFFERENCE  FOUND. 


COMPUTE  THE  MAXIMUM  DIFFERENCE  OF  THE  CURRENT 
CONFIGURATION. 

NBW  =  0 

DO  200  N  =  1,  NUMEL 
DO  100  J  =  1,  7 
J1  =  IX(J,  N) 

KK  =  J  +  1 
DO  60  K  =  KK,  8 
K1  =  IX(K,  N) 

KD  =  IABS(J1  -  Kl) 

NBW  =  MAX0(NBW,  KD) 

CONTINUE 
CONTINUE 
CONTINUE 

PRINT  205,  NBW 

FORMAT  (/  'INITIAL  MAXIMUM  DIFFERENCE  =  ',  15) 

DETERMINE  THE  SMALLEST  MAXIMUM  DIFFERENCE 
POSSIBLE. 

NBWMIN  =  NBW 

ITERATE  TO  A  SOLUTION. 

NOITER  =  MINC(NUMNP,  20) 

DO  220  N  =  1,  NUMNP 
NODE(N)  =  N 
220  CONTINUE 

DO  1100  L  =  1,  NOITER 

INITIALIZE. 

DO  300  N  =  1,  NUMNP 

NNS(N)  =  99999 

NOS{N)  =  0 

300  CONTINUE 
C 

K  =  1 

NNS(l)  =  L 

NOS(L)  =  1 

NDIFF  =  0 


C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 


60 

100 

200 


205 


186 


ooo  ooonooooo  oooo  oo  non 


PROCESS  THE  ADJACENCY  TABLE. 

DO  800  N  =  1,  NUMNP 

ISUB  =  NNS(N) 

NUM  =  JZ(ISUB,  9) 

DO  700  M  =  1,  NUM 

NBR  =  JZ(ISUB,  M) 

IF  (NOS(NBR))  400,  400,  700 
400  K  =  K  +  1 

NOS (NBR)  =  K 
NNS(K)  =  NBR 

CHECK  IF  MINIMUM  MAXIMUM  DIFFERENCE  HAS  BEEN 
EXCEEDED. 

KD  =  IABS(N  -  K) 

NDIFF  =  MAX0(NDIFF,  KD) 

IF  (KD  -  NBWMIN)  500,  500,  1100 

CHECK  IF  ITERATION  IS  FINISHED. 

500  IF  (K  -  NUMNP)  700,  900,  900 

700  CONTINUE 

800  CONTINUE 

DETERMINE  THE  MAXIMUM  DIFFERENCE  OF  THE  CURRENT 
CONFIGURATION. 

900  NBW  =  NDIFF 


DO 

920  N  =  1, 

NUMEL 

DO 

910  J  =  1, 

7 

J1 

=  IX(J,  N) 

J1 

=  NOS(Jl) 

KK 

=  J  +  1 

DO 

905  K  =  KK, 

8 

K1 

=  IX(K,  N) 

K1 

=  NOS(Kl) 

KD 

=  IABS(J1  - 

Kl) 

NBW  =  MAX0(NBW 

,  KD) 

905  CONTINUE 
910  CONTINUE 
920  CONTINUE 

CHECK  FOR  THE  BEST  CONFIGURATION. 

IF  (NBWMIN  -  NBW)  1100,  1100,  950 
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950  NBWMIN  =  NBW 

DO  1000  N  =  1,  NUMNP 
NODE(N)  =  NOS(N) 

1000  CONTINUE 
C 

PRINT  1010,  L,  NBWMIN 

1010  FORMAT  ('ITERATION  =  15,  lOX,  'BEST  MAXIMUM' 

&  '  DIFFERENCE  =  ' ,  15) 

C 

1100  CONTINUE 

SHUFFLE  THE  DATA. 

DO  1200  N  =  1,  NUMNP 
N1  =  NODE(N) 

NNS(Nl)  =  N 
1200  CONTINUE 

FIX  THE  NODE  DATA. 

REWIND  8 

DO  1400  N  =  1,  NUMNP 
N1  =  NNS(N) 

WRITE  (8)  X(N1),  Y(N1),  Z(N1),  IBC(Nl),  BV(Nl) 
1400  CONTINUE 
REWIND  8 

DO  1500  N  =  1,  NUMNP 

READ  (8)  X(N),  Y(N),  Z(N),  IBC(N),  BV(N) 

1500  CONTINUE 

FIX  THE  ELEMENT  DATA. 

DO  1700  N  =  1,  NUMEL 
DO  1600  I  =  1,  8 

11  =  IX(I,  N) 

12  =  NODE(Il) 

IX(I,  N)  =  12 

1600  CONTINUE 
1700  CONTINUE 
C 

RETURN 

END 


The  most  difficult  aspect  of  this  subroutine  is  keeping 
up  with  the  new  node  numbers  (the  two  arrays  NOS  and  NNS 
accommodate  this) .  Finally,  any  number  of  tries  can  be  made 
as  the  starting  point  to  renumber  the  grid  up  to  the  number 
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of  node  points.  This  subroutine  currently  has  the  number  of 
iterations  as  20. 


Subroutine  BC 

Subroutine  BC  applies  boundary  conditions,  and  its 
listing  is  now  given: 


C 

C 

C 

C 

C 

C 


C 

C 

C 

C 


C 

C 

C 

C 

C 

c 


c 


SUBROUTINE  BC(IBLN0,  IIB,  JIB,  KIB,  I2B,  J2B,  K2B, 
&  NBC,  BCVAL) 


THIS  SUBROUTINE  APPLIES  BOUNDARY  CONDITIONS  TO  THE 
FACES  OF  BLOCKS. 


PARAMETER  (ND  =  1000) 

COMMON  /OUTPUT/  X(ND),  Y(ND),  Z(ND),  IBC(ND),  BV(ND), 

&  IX(9,  ND) ,  THETA(ND),  PHI(ND),  PSI (ND) ,  NUMNP,  NUMEL, 
&  NUMMAT,  NODVEL 

DIMENSION  IPICK(2,  2,  2),  IFACE(4,  6),  JNODE(4) 

DATA  IPICK  /I,  2,  4,  3,  5,  6,  8,  7/ 

DATA  IFACE  /I,  5,  8,  4,  2,  3,  7,  6,  1,  2,  6,  5,  3,  4, 

&  8,  7,  4,  3,  2,  1,  5,  6,  7,  8/ 

GET  SIZE  OF  BLOCK  AND  NUMBER  OF  ELEMENTS  JUST 
BEFORE  IT  BEGINS. 

READ  (7,  REC=IBLNO)  12,  J2 ,  K2 ,  NEL 

IROW  =12-1 

IPLANE  =  (J2  -  1)  *  IROW 

APPLY  VARIOUS  BOUNDARY  CONDITIONS,  EXCEPT  NBC  = 

-2  . 

IF  (NBC.EQ.-2)  GO  TO  3050 

DO  3000  K  =  KIB,  K2B 

IF  (K.LT.K2)  THEN 
KK  =  1 
ELSE 

KK  =  2 
ENDIF 


C 


DO  2900  J  =  JIB,  J2B 
IF  (J.LT.J2)  THEN 
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JJ  =  1 
ELSE 

JJ  =  2 
ENDIF 

DO  2800  I  =  IIB,  I2B 

IF  (I.LT.I2)  THEN 
II  =  1 
ELSE 

II  =  2 
ENDIF 

GET  NODE  NUMBER. 

MM  =  (K  -  KK)  *  IPLANE  +  (J  -  JJ)  *  IROW  +  I  - 
&  +  NEL 

IPOS  =  IPICK(II,  JJ,  KK) 

NN  =  IX (IPOS,  MM) 

DO  VARIOUS  BOUNDARY  CONDITIONS. 

IF  ( (NBC.EQ.l) .OR. (NBC.EQ.-l) )  THEN 
IBC(NN)  =  NBC 
BV(NN)  =  BCVAL 

ELSE  IF  ( (NBC.EQ.2) .OR. (NBC.EQ. 12) )  THEN 
IBC(NN)  =  NBC 
BV(NN)  =  Z(NN) 

ENDIF 

2800  CONTINUE 
2900  CONTINUE 
3000  CONTINUE 

GO  TO  5000 

DO  SPECIFIED  DISCHARGE  VELOCITY  FOR  A  FACE 

3050  IF  ( (IIB.EQ. 1) . AND. (I2B.EQ. 1) )  THEN 
JF  =  1 
LI  =  JIB 
L2  =  J2B  -  1 
Ml  =  KIB 
M2  =  K2B  -  1 
LFACT  =12-1 
MFACT  =  LFACT  *  (J2  -  1) 

MST  =  NEL  +  1 

ELSE  IF  ( (IIB.EQ. 12) .AND. (I2B.EQ. 12) )  THEN 
JF  =  2 
LI  =  JIB 


II  +  1 
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1 


L2  =  J2B  - 


Ml  =  KIB 

M2  =  K2B  -  1 

LFACT  =12-1 

MFACT  =  LFACT  *  (J2  -  1) 

MST  =  NEL  +12-1 

ELSE  IF  ( (JIB.EQ. 1) .AND. (J2B.EQ. 1) )  THEN 
JF  =  3 
LI  =  IIB 
L2  =  I2B  -  1 
Ml  =  KIB 
M2  =  K2B  -  1 
LFACT  =  1 

MFACT  =  (12  -  1)  *  (J2  -  1) 

MST  =  NEL  +  1 

ELSE  IF  ( (JIB.EQ. J2) .AND. (J2B.EQ.J2) )  THEN 
JF  =  4 
LI  =  IIB 
L2  =  I2B  -  1 
Ml  =  KIB 
M2  =  K2B  -  1 
LFACT  =  1 

MFACT  =  (12  -  1)  *  (J2  -  1) 

MST  =  NEL  +  (12  -  1)  *  (J2  -  2)  +  1 

ELSE  IF  ( (KIB.EQ. 1) .AND. (K2B.EQ. 1) )  THEN 
JF  =  5 
LI  =  IIB 
L2  =  I2B  -  1 
Ml  =  JIB 
M2  =  J2B  -  1 
LFACT  =  1 
MFACT  =12-1 
MST  =  NEL  +  1 

ELSE  IF  ( (KIB.EQ. K2) .AND. (K2B.EQ.K2) )  THEN 
JF  =  6 
LI  =  IIB 
L2  =  I2B  -  1 
Ml  =  JIB 
M2  =  J2B  -  1 
LFACT  =  1 
MFACT  =12-1 

MST  =  NEL  +  (12  -  1)  *  (J2  -  1)  *  (K2  -  2)  +  1 

ENDIF 

DO  3500  M  =  Ml,  M2 

DO  3400  L  =  LI,  L2 

GET  NODE  NUMBERS. 

1)  *  LFACT  +  (M  - 
=  1,  4 
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MM  =  (L  - 
DO  3100  I 


1)  *  MFACT  +  MST 


JPOS  =  IFACE(I,  JF) 

JJ  =  IX (JPOS,  MM) 
IBC(JJ)  =  -1 
JNODE(I)  =  JJ 
3100  CONTINUE 
C 

C  SAVE  RESULTS . 

C 

NODVEL  =  NODVEL  +  1 
WRITE  (9)  JNODE,  BCVAL 
C 

3400  CONTINUE 
C 

3500  CONTINUE 
C 

5000  RETURN 
END 


This  subroutines  is  roughly  divided  into  two  parts. 

The  first  part  applies  boundary  conditions  1,  2,  12,  or  -1 
to  a  group  on  nodes,  typically  a  sheet  on  the  surface  of  the 
block.  The  last  section  determines  sets  of  four  nodes 
required  for  discharge  velocity  data  and  writes  them  to  a 
temporary  file.  Since  the  elements  are  not  shuffled  in 
bandwidth  minimization,  the  proper  nodes  can  be  obtained  by 
looking  up  the  correct  node  of  the  correct  element  in  the 
connectivity  array.  This  process  is  greatly  sim  /..ified  by 
the  arrays  IPICK  and  IFACE. 


Subroutine  COMBIN 

Subroutine  COMBIN  reads  in  any  nr.nber  of  EAGLE  grid 
files  and  combines  them  into  a  single  FEM  grid.  Its  listing 
follows . 
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SUBROUTINE  COMBIN (NOFTLE ,  IBL) 

C 

C 

C  THIS  SUBROUTINE  TAKES  TWO  OR  MORE  FILES  CONTAINING 

C  EAGLE  GRIDS  AND  COMBINES  THEM  INTO  ONE  GRID. 

C 

C 

PARAMETER  (ND  =  1000) 

COMMON  /OUTPUT/  X(ND),  Y(ND),  Z (ND) ,  IBC(ND),  BV(ND), 

&  IX(9,  ND) ,  THETA (ND) ,  PHI(ND),  PSI(ND),  NUMNP,  NUMEL, 
&  NUMMAT,  NODVEL 

COMMON  /SCRAT/  NOD(ND),  IREPL(ND) ,  BB(ND) 

CHARACTER  ANAME  *  20 

DIMENSION  I2S(20),  J2S(20),  K2S(20) 

C 

IBL  =  0 
IZERO  =  0 

LOOP  OVER  THE  NUMBER  OF  FILES. 

DO  1000  M  =  1,  NOFILE 

READ  EAGLE  GRID  FILE. 

10  PRINT  15,  M 

15  FORMAT (/  '  FILE  NAME  NO.',  13,  '  ?') 

20  READ  (*,  25)  ANAME 
25  FORMAT (A) 

OPEN  (2,  FILE=ANAME,  STATUS= ' OLD ' ,  IOSTAT=ISTAT) 

IF  (ISTAT.NE.O)  THEN 
PRINT*,  '  ' 

NN  =  NONBLK( ANAME) 

PRINT*,  'INCORRECT  FILE  NAME  -  ',  ANAME ( 1  :  NN) , 

&  '  -  TRY  AGAIN. ' 

GO  TO  20 
ENDIF 

GET  THE  NUMBER  OF  BLOCKS. 

READ  (2,  *)  NUMBLK 
IF  (NUMBLK. GT. 20)  THEN 

PRINT*,  'NO  MORE  THAN  20  BLOCKS  MAY  BE  IN', 

&  '  ONE  FILE. ' 

STOP 
ENDIF 

GET  THE  GRID  DIMENSIONS  FOR  ALL  BLOCKS. 

DO  30  L  =  1,  NUMBLK 
READ  (2,  *)  I2S(L),  J2S(L),  K2S(L) 

30  CONTINUE 
C 


193 


c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 

c 


c 


c 


LOOP  OVER  THE  NUMBER  OF  BLOCKS. 

DO  900  L  =  1,  NUMBLK 

GET  THE  SIZE  OF  THE  BLOCK. 

12  =  I2S(L) 

J2  =  J2S(L) 

K2  =  K2S(L) 

IJ2  =  12  *  J2 

OUTPUT  BLOCK  DATA  TO  TEMPORARY  FILE. 

IBL  =  IBL  +  1 

WRITE  (7,  REC=IBL)  12,  J2,  K2 ,  NUMEL 

READ  THE  (X,  Y,  Z)  DATA  POINTS. 

NUM  =  12  *  J2  *  K2 
NDl  =  NUMNP  +  1 
NDTOT  =  NUMNP  +  NUM 

READ  (2,  *)  (X(I),  I  =  NDl,  NDTOT),  (Y(I),  I  =  NDl, 

&  NDTOT),  (Z(I),  I  =  NDl,  NDTOT) 

COMPUTE  ELEMENT  DATA. 

NN  =  NUMNP 
MM  =  NUMEL 

GET  MATERIAL  TYPE  NUMBER,  THETA,  PHI,  AND  PSI. 
PRINT  42,  L 

42  FORMAT(/  '  FOR  BLOCK',  13,  '  MATERIAL  TYPE  NUMBER, 

&  THETA,  PHI,',  '  AND  PSI?') 

READ*,  MAT,  TH,  PH,  PS 
NUMMAT  =  MAXO(NUMMAT,  MAT) 


DO  70 

K  = 

1, 

K2  -  1 

DO  60 

J  = 

1, 

J2  -  1 

DO  50 

I  = 

1, 

12-1 

NN  = 

NN  + 

1 

MM  = 

MM  + 

1 

IX(1, 

MM) 

= 

NN 

IX(2, 

MM) 

= 

IX(1,  MM) 

+  1 

IX(3, 

MM) 

= 

IX(2,  MM) 

+  12 

IX(4, 

MM) 

= 

IX(3,  MM) 

-  1 

DO  45  II  =  1,  4 
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IX(II+4,  MM)  =  IX(II,  MM)  +  IJ2 
45  CONTINUE 
C 

IX (9,  MM)  =  MAT 
THETA (MM)  =  TH 
PHI  (MM)  =  PH 
PSI(MM)  =  PS 
C 

50  CONTINUE 
C 

NN  =  NN  +  1 
C 

60  CONTINUE 
C 

NN  =  NN  +  12 
C 

70  CONTINUE 
C 

NUMEL  =  MM 

INITIATE  THE  NODE  ARRAY. 

DO  220  I  =  .1,  NDTOT 
NOD (I)  =  I 
IREPL(I)  =  0 
220  CONTINUE 

REMOVE  DUPLICATE  NODES. 

IF  (L.EQ.l)  THEN 
N1  =  1 
ELSE 

N1  =  NDl 
ENDIF 
C 

DO  300  N  =  Nl,  NDTOT 
C 

XCK  =  X(N) 

YCK  =  Y(N) 

ZCK  =  Z(N) 

IF  (MOD(N,  1000). EQ.O)  PRINT*,  'REMOVING 
&  ' NODES ,  N  = ' ,  N ,  '  NUMNP  = ' ,  NUMNP 

C 

DO  240  1=1,  N-1 
IF  (IREPL(I) .EQ. 1)  GO  TO  240 
CK  =  ABS(X(I)  -  XCK)  +  ABS(Y(I)  -  YCK)  + 
&  ZCK) 

IF  (CK  -  l.E-2)  250,  250,  240 
240  CONTINUE 
GO  TO  300 
C 

250  IREPL(N)  =  1 


DUPLICATE  ' 


ABS(Z(I)  - 
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NOD(N)  =  NOD (I) 

IF  (N.EQ.NDTOT)  GO  TO  300 
C 

260  NPl  =  N  +  1 

DO  270  I  =  NPl,  NDTOT 
NOD(I)  =  NOD(I)  -  1 
270  CONTINUE 
C 

300  CONTINUE 

UPDATE  THE  NODE  DATA. 

LL  =  NUMNP 

DO  320  I  =  NDl,  NDTOT 
IF  (IREPL(I) .EQ. 1)  GO  TO  320 
LL  =  LL  +  1 
X(LL)  =  X(I) 

Y(LL)  =  Y(I) 

Z(LL)  =  Z(I) 

IBC(LL)  =  IBC(I) 

BV(LL)  =  BV(I) 

320  CONTINUE 

NUMNP  =  LL 

UPDATE  THE  ELEMENT  DATA. 

DO  340  1=1,  NUMEL 
DO  330  J  =  1,  8 
INOD  =  IX(J,  I) 

IX(J,  I)  =  NOD(INOD) 

330  CONTINUE 
340  CONTINUE 
C 

900  CONTINUE 
C 

CLOSE  (2) 

C 

1000  CONTINUE 
C 

RETURN 

END 


After  the  EAGLE  data  is  read  into  memory,  the  element 
connectivity  data  is  then  created.  Next,  duplicate  nodes 
are  removed  by  considering  distance  squared  between  two 
given  node  points.  When  nodes  have  been  removed,  the 


remaining  node  points  are  renumbered  so  they  are 
consecutively  numbered.  With  this  done,  the  element  connec¬ 
tivity  must  then  be  modified  to  reflect  the  new  node 
numbers . 


Subroutine  DONE 

Subroutine  DONE  closes  all  files  that  are  still  open, 
and  its  listing  is  given. 


SUBROUTINE  DONE 
C 
C 

C  THIS  SUBROUTINE  WRAPS  THINGS  UP. 

C 

C 

CLOSE  (7) 

CLOSE  (8) 

CLOSE  (9) 

C 

RETURN 

END 


Subroutine  OUTFEM 

Subroutine  OUTFEM  writes  the  REM  grid  in  3-D  FEM 
seepage/groundwater  format.  Its  listing  is  as  follows: 


SUBROUTINE  OUTFEM (IBL) 

C 

C 

C  THIS  SUBROUTINE  OUTPUTS  A  FILE  FOR  THE  3-D  SEEPAGE 

C  GROUNDWATER  FEM  PROGRAM. 

C 

C 

PARAMETER  (ND  =  1000) 

COMMON  /OUTPUT/  X(ND),  Y(ND),  Z(ND),  IBC(ND),  BV(ND), 

&  IX(9,  ND) ,  THETA(ND),  PHI(ND),  PSI(ND),  NUMNP,  NUMEL, 
&  NUMMAT,  NODVEL 
CHARACTER  TITLE  *  80 
DIMENSION  JNODE(4) 

C 

C  DO  TITLE. 
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c 


PRINT*,  '  ' 

PRINT*,  'TITLE  CARD?' 

READ  100,  TITLE 
100  FORMAT(A80) 

WRITE  (2,  100)  TITLE 
C 

C  DO  NUMBER  OF  NODE  POINTS,  NUMBER  OF  ELEMENTS, 

C  NUMBER  OF  DIFFERENT  MATERIALS,  NUMBER  OF  DISCHARGE 

C  VELOCITY  CARDS,  AND  DATUM. 

C 

PRINT* ,  '  ' 

PRINT* ,  ' DATUM? ' 

READ  (*,  *)  DATUM 

WRITE  (2,  200)  NUMNP,  NUMEL,  NUMMAT,  NODVEL,  DATUM 
200  FORMAT(4I5,  F10.2) 

C 

C  ADJUST  BOUNDARY  CONDITION  DATA  FOR  IBC  =  2  AND 

C  IBC  =  12. 

C 

DO  205  1=1,  NUMNP 

IF  ( (IBC(I) .EQ.2) .OR. (IBC(I) .EQ. 12) )  BV ( I )  =  BV(I)  - 
&  DATUM 
205  CONTINUE 

DO  MATERIAL  DATA. 

DO  250  1=1,  NUMMAT 
PRINT  210,  I 

210  FORMAT (/  '  Kl,  K2 ,  AND  K3  FOR  MATERIAL  TYPE',  13, 

&  '  ?  '  ) 

READ  (*,  *)  FKl,  FK2,  FK3 
WRITE  (2,  220)  I,  FKl,  FK2 ,  FK3 
220  FORMAT(I5,  3E10.3) 

250  CONTINUE 

DO  NODES. 

IF  (NUMNP. EQ.O)  GO  TO  1000 
DO  400  I  =  1,  NUMNP 

WRITE  (2,  300)  I,  IBC(I),  X(I),  Y(I),  Z(I),  BV(I) 

300  FORMAT(2I5,  4F10.2) 

400  CONTINUE 

DO  ELEMENTS. 

DO  600  J  =  1,  NUMEL 

WRITE  (2,  500)  J,  (IX(I,  J),  1=1,  9),  THETA(J) 

&  PHI (J) ,  PSI (J) 

500  FORMAT(10I5,  3F10.2) 

600  CONTINUE 


C 
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C  DO  DISCHARGE  VELOCITY  DATA. 

C 

IF  (NODVEL.EQ. 0)  GO  TO  900 
REWIND  9 

DO  800  J  =  1,  NODVEL 
READ  (9)  JNODE,  BCVAL 
WRITE  (2,  700)  JNODE,  BCVAL 
700  FORMAT(4I5,  F10.2) 

800  CONTINUE 
C 

C  DO  GRID  DIMENSION  DATA. 

C 

900  WRITE  (3,  200)  IBL 
DO  910  1=1,  IBL 

READ  (7,  REC  =  I^  12,  J2,  K2 ,  NEL 
WRITE  (3,  200)  12,  J2 ,  K2 
910  CONTINUE 
C 

1000  RETURN 
END 


Additional  information  is  obtained  as  needed  as  the 
output  process  is  done.  This  includes  a  title  line,  the 
datum,  and  material  data.  Notice  that  the  head  values  are 
corrected  by  the  datum  if  boundary  condition  2  or  12  is 
specified . 


Function  NONBLK 

This  function  counts  the  number  of  characters  in  a  file 
name  so  the  extensions,  .egl,  .eg2,  .dim,  .sol,  etc.  can  be 
properly  attached.  Its  listing  is  now  given: 


FUNCTION  NONBLK ( FI LENM) 

C 

C 

C  THIS  FUNCTION  COMPUTES  THE  NUMBER  OF  NON-BLANK 

C  CHARACTERS  IN  THE  FILE  NAME  FI LENM,  IT  WORKS 

C  BY  LOOKING  FOR  THE  FIRST  BLANK  CHARACTER. 

C 

C 

CHARACTER  FILENM  *  (*) 


C 
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N  =  0 

DO  100  1=1,  100 

IF  (FILENM(I  :  I).EQ.'  •)  GO  TO  200 

N  =  N  +  1 
100  CONTINUE 
C 

200  NONBLK  =  N 
RETURN 
END 
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APPENDIX  C 

CONVERSION  TO  FAST  FORMAT 
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This  appendix  contains  the  subroutine  to  convert  the 


output  from  the  unstructured  finite  element  grid  to  the 
multiblocked  structured  finite  volume  grid  format  used  by 
FAST.  Although  FAST  was  written  primarily  to  display  compu 
tational  fluid  dynamics  (CFD)  data,  it  can  also  be  used  to 
present  other  phenomena  as  well.  The  listing  of  the  subrou 
tine  is  now  given. 


SUBROUTINE  FASTOT ( DATUM) 

C 

C 

C  THIS  SUBROUTINE  CALCULATES  THE  POSTPROCESSOR  FILE 

C  FOR  FAST . 

C 

C 

PARAMETER  (NDX  =  25000) 

PARAMETER  (ND  =  1250) 

PARAMETER  (ND2  =  ND  *  2) 

COMMON  /  GRIDl  /  NUMNP,  NUMEL,  NUMMAT,  X(NDX),  Y(NDX) 
&  Z(NDX),  FX(NDX),  NBC(NDX),  FLOW(NDX),  HLAST(NDX) 
COMMON  /  GRID2  /  XK1(12),  XK2(12),  XK3(12),  NP(9, 

&  NDX),  THETA (NDX),  PHI (NDX) ,  PSI(NDX) 

COMMON  /  FREE  /  FRX(NDX),  FRY(NDX),  FRZ(NDX), 

&  COUNT (NDX) 

COMMON  /  BANARG  /  MBAND,  NUMBLK,  R(NDX),  C(ND2,  ND) 
DIMENSION  VX (NDX) ,  VY(NDX),  VZ(NDX),  Q(NDX,  8) 
EQUIVALENCE  (FRX,  VX) ,  (FRY,  VY) ,  (FRZ,  VZ),  (C,  Q) 
DIMENSION  IPICK(2,  2,  2) 

DATA  IPICK  /I,  2,  4,  3,  5,  6,  8,  7/ 

C 

C  COMPUTE  MAXIMUM  AND  MINIMUM  HEAD. 

C 

HMIN  =  1.E30 

HMAX  =  -  HMIN 

DO  100  1=1,  NUMNP 

HMIN  =  AMIN1(R(I)  -  DATUM,  HMIN) 

HMAX  =  AMAX1(R(I)  -  DATUM,  HMAX) 

100  CONTINUE 

RDH  =  1.  /  (HMAX  -  HMIN) 

C 

C  SET  AERODYNAMIC  DATA. 

C 

FSMACH  =  HMIN 
ALPHA  =  HMAX 
RE  =  1. 
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TIME  =  0. 


WRITE  NUMBER  OF  BLOCKS. 

READ  (3,  *)  IBL 
WRITE  (12,  130)  IBL 
WRITE  (13,  130)  IBL 
130  FORMAT(3I5) 

WRITE  GRID  DIMENSION  DATA. 

DO  150  N  =  1,  IBL 
READ  (3,  *)  12,  J2,  K2 
WRITE  (12,  130)  12,  J2,  K2 
WRITE  (13,  130)  12,  J2 ,  K2 
150  CONTINUE 


300 


C 

C 

C 

C 

C 

C 


WRITE  SOLUTION  FOR  ALL  GRIDS. 

REWIND  3 

READ  (3,  *)  IBL 

NEL  =  0 

DO  1000  N  =  1,  IBL 

READ  (3,  *)  12,  J2,  K2 

IROW  =12-1 

IPLANE  =  (J2  -  1)  *  IROW 

WRITE  PRELIMINARY  DATA. 

WRITE  (12,  300)  FSMACH,  ALPHA,  RE,  TIME 
FORMAT(6(lX,  Ell, 5)) 

ACCUMULATE  Q  VECTOR  WITH  RESULTS  IN  UK  ORDER. 
INOD  =  0 

DO  800  K  =  1,  K2 
KK  =  1 

IF  (K.EQ.K2)  KK  =  2 
DO  700  J  =  1,  J2 
JJ  =  1 

IF  (J.EQ.J2)  JJ  =  2 
DO  600  I  =  1,  12 
II  =  1 

IF  (I.EQ.I2)  II  =  2 
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GET  NODE  NUMBER. 

MM  =  (K  -  KK)  *  IPLANE  +  (J  -  JJ)  *  IROW  +  I  -  II  +  1 
&  +  NEL 

IPOS  =  IPICK(II,  JJ,  KK) 

NN  =  NP(IPOS,  MM) 

FOR  NODES  ABOVE  THE  FREE  SURFACE  SET  HEAD  EQUAL  TO 
ELEVATION. 

H  =  R(NN) 

IF  (H.LT.Z(NN))  H  =  Z(NN) 

NORMALIZE  HEAD  BETWEEN  1  AND  2. 

HN  =  (H  -  DATUM  -  HMIN)  *  RDH  +  1. 

COMPUTE  Q. 

INOD  =  INOD  +  1 
Q(INOD,  1)  =  HN 

Q(INOD,  2)  =  VX(NN)  *  HN 

Q(INOD,  3)  =  VY(NN)  *  HN 

Q(INOD,  4)  =  VZ(NN)  *  HN 

Q(INOD,  5)  =  H  -  DATUM 

USE  THE  LAST  THREE  POSITIONS  FOR  Q  FOR  THE  NEW 
(X,  Y,  Z). 

Q(INOD,  6)  =  X(NN) 

Q(INOD,  7)  =  Y(NN) 

Q(INOD,  8)  =  Z(NN) 

600  CONTINUE 

700  CONTINUE 

800  CONTINUE 


OUTPUT  Q. 


WRITE  (12, 
WRITE  (13, 
C 

NEL  =  (12 
C 

1000  CONTINUE 
C 

CLOSE  (7) 

C 


300)  ( (Q(I,  J) ,  1=1, 

300)  ((Q(I,  J),  1=1, 

-  1)  *  (J2  -  1)  *  (K2  - 


INOD)  ,  J 
INOD)  ,  J 

1)  +  NEL 


1,  5) 

6,  8) 
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RETURN 


END 

Because  in  an  unconfined  flow  problem  the  original  grid 
is  modified,  both  a  file  containing  the  new  grid  and  a  file 
containing  the  results  are  written.  C  programs  were  written 
to  take  the  ASCII  output  files  from  this  subroutine  and 
convert  them  to  binary  files  (FORTRAN  binary  WRITE ' s  will 
not  work) .  The  C  programs  could,  of  course,  be  converted  to 
functions  and  called  from  FASTOT  as  well.  Because  a  totally 
different  flow  is  being  modeled,  the  preliminary  data  such 
as  free  stream  Mach  Number,  angle  of  attack,  and  Reynold's 
Number  are  used  for  other  things.  Mach  Number  has  minimum 
head  stored  in  its  place,  angle  of  attack  has  maximum  head, 
and  Reynold's  Number  is  set  to  1.  The  unconfined  flow  part 
of  the  seepage/groundwater  solver  is  a  solution  to  a  system 
of  nonlinear  equations  rather  than  using  time  dependent 
equations  and  allowing  the  computation  to  converge  to  a 
steady-state  condition.  Therefore,  integration  time  is  not 
known,  and  it  is  set  to  zero. 

The  remaining  data  to  be  written  to  the  output  file  are 
components  of  the  Q'  vector 
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where 

p  =  density 

u  =  X  component  of  velocity 
v  -  y  component  of  velocity 
w  =  z  component  of  velocity 
U  =  energy 

Again,  because  of  the  difference  in  applications,  total  head 
instead  of  density  and  discharge  velocity  instead  of 
velocity  are  used.  In  fact,  head  is  normalized  between  one 
and  two  as  follows; 

h  =  —  ~  +  1 

An  ax  -  Anxn 

where 

h^  =  normalized  head 
^min  ~  minimum  head 
=  maximum  head 

This  normalization  prevented  division  by  zero  by  FAST  in 
computing  the  components  of  velocity 
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The  momentum  type  components  are  then  written  as  h^lv), 
where  {v}  is  the  discharge  velocity  vector.  Energy  is  not 
used,  so  the  original  values  of  head  are  output  instead. 

The  array  IPICK  makes  it  significantly  easier  to  choose 
one  of  the  eight  nodes  of  a  given  element  on  which  to  write 
information  (see  Chapter  V  under  the  heading,  "Scientific 
Visualization,"  and  the  subheading,  "Conversion  from  Finite 
Element  to  Finite  Volume  Format,"  for  more  details). 

Note  also  that  all  values  of  the  Q'  vector  are  first 
computed  and  then  each  component  is  written  to  the  output 
file.  This  is  because  FAST  format  (originally  PL0T3D) 
requires  all  values  of  a  given  component  of  the  Q  vector  to 
be  written  before  another  component  is  processed. 
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APPENDIX  D 
GRID  GENERATION  DATA 


This  appendix  gives  the  data  required  to  generate  the 
elliptic  grid  for  the  tvvo-well-in-an-aquifer  problem  given 
in  Chapter  VI.  These  include  two  O  blocks  generated  by 
EAGLE,  two  plugs  generated  by  EAGLE,  and  the  data  to  combine 
and  prepare  these  data  into  FEM  seepage/groundwater  format. 


Surface  Generation  Data 

The  data  for  the  surface  generation  portion  of  EAGLE 
for  the  first  well  O  grid  is  now  given. 


$  'point',  point=l,  r=175 . , -200 . , 0 .  $ 

$  'point',  point=2,  r=175 . , 300 . , 0 .  $ 

$  'point',  point=3,  r=-250 . , 300 . ,  0 .  $ 

$  'point',  point=4,  r=-250 . , -200 . , 0 .  $ 

$  'point',  point=5,  r=175 . , -200 . , 120 .  $ 

$  'point',  point=6,  r=175 . , 300 .  ,  120 .  $ 

$  'point',  point=7 ,  r=-250. , 300. , 120.  $ 

$  'point',  point=8,  r=-250 . , -200 . ,  120 .  $ 

$  'setnum',  segment=l,  points=13  $ 

$  'setnum',  segment=2 ,  points=6  $ 

$  'setnum',  segment=3 ,  points=25  $ 

$  'setnum',  segment=4 ,  points=9  $ 

$  'setnum',  segment=5,  points=17  $ 

$  'setval',  number=l,  value=.002  $ 

$  'setval',  number=2,  value=.04  $ 

$  'setval',  number=3,  value=.02  $ 

$  'conicur',  points=-2,  type= ' circle ' ,  radius=l., 
angle=59 . 7 , -48 . 8 ,  coreout=5  $ 

$  'conicur',  points=-2,  type= ' circle ' ,  radius=l., 
angle=129 . 8 , 59 . 7 ,  coreout=6  $ 

$  'conicur',  points=-2,  type= ' circle ' ,  radius=l., 
angle=218 . 7 , 129 . 8 ,  coreout=7  $ 

$  'conicur',  points=-2,  type= ' circle ' ,  radius=l., 
angle=-48 . 8 , -14 1 . 3 ,  coreout=8  $ 

$  'trans',  corein=5,  origin=0 . , 0 . , 80 . ,  coreout=49  $ 

$  'trans',  corein=5,  origin=0 . , 0 . , 120 . ,  coreout=45  $ 
$  'getend',  corein=49,  point= ' first ' ,  end=' first'  $ 

$  'getend',  corein--45,  point- '  first '  ,  end='last'  $ 

$  'line',  points=-4,  space=-2,  coreout=13  $ 

$  'getend',  corein-49,  point= ' f irst '  ,  end=' first'  $ 

$  'getend',  corein=5,  point= ' first ' ,  end='last'  $ 

$  'line',  points=-5,  space=-3  $ 

$  'switch',  reorder= ' reversel '  $ 

$  'insert',  coreln=13,  coreout=13  $ 
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$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 


'line',  points=-3,  rl=2,  r2=6,  coreout=17  $ 

'getend',  corein=5,  point= ' f irst ' ,  end=' first'  $ 

'line',  points=-l,  r2=2,  space=-l,  coreout=2  $ 

'trans',  corein=2,  origin=0 . , 0 . , 120 . ,  coreout=26  $ 

'  edgecur ' ,  edge= ' lowerl ' ,  corein=13  $ 

'edgecur',  edge= ' upperl ' ,  corein=17  $ 

'edgecur',  edge= ' lower2 * ,  corein=2  $ 

'edgecur',  edge= ' upper2 ' ,  corein=26  $ 

'transur',  coreout=51  $ 

'current',  corein=51,  coreout=52  $ 

'bouncur',  corein=13  $ 

'rotate',  angpts=-2,  angle=0 . , -108 . 5 ,  axcos=0 . , 0 . , 1 .  $ 
'switch',  reorder= ' switch ' ,  coreout=41  $ 

'bouncur',  corein=13  $ 

'rotate',  angpts=-2 ,  angle=-108 . 5 , -201 . ,  axcos=0 . , 0 . , 1 .  $ 
'switch',  reorder= ' switch ' ,  coreout=42  $ 

'bouncur',  corein=13  $ 

'rotate',  angpts=-2 ,  angle=-201 . , -289 . 9 ,  axcos=0 . , 0 . , 1 .  $ 
'switch',  reorder= ' switch ' ,  coreout=43  $ 

'bouncur',  corein=13  $ 

'rotate',  angpts=-2 ,  angle=70 . 1 , 0 . ,  axcos=0 . , 0 . , 1 .  $ 
'switch',  reorder= ' switch ' ,  coreout=44  $ 

'current',  corein=41  $ 

'insert',  corein=42  $ 

'insert',  corein=43  $ 


insert 

'  ,  corein=44 

,  coreout=53 

i  $ 

1  ine '  , 

points= 

-3, 

rl=l. 

r2  =  5. 

coreout=20 

$ 

line '  , 

points= 

-3, 

rl  =  2. 

r2=6. 

coreout=17 

$ 

line '  , 

points= 

-3, 

rl=3, 

r2=7, 

coreout=18 

$ 

1  ine '  , 

points= 

-3, 

rl=4, 

r2=8, 

coreout=19 

$ 

line '  , 

points= 

-2, 

rl=2 , 

r2  =  l. 

coreout=9 

$ 

line '  , 

points^ 

-2, 

rl  =  3 , 

r2=2. 

coreout=10 

$ 

line '  , 

points= 

-2, 

rl=4. 

r2=3. 

coreout=ll 

$ 

line '  , 

points= 

-2, 

rl=l. 

r2=4 , 

coreout=12 

$ 

line '  , 

points= 

-2, 

rl=6 , 

r2=5. 

coreout=22 

$ 

line '  , 

points= 

-2, 

II 

r2=6. 

coreout=2  3 

$ 

line '  , 

points^ 

-2, 

rl=8. 

r2=7. 

coreout=24 

$ 

line '  , 

points= 

-2, 

rl=5. 

r2=8. 

coreout=2 1 

$ 

'edgecur',  edge= ' lowerl ' ,  corein=17  $ 
'edgecur',  edge= ' upperl ' ,  corein=20  $ 
'edgecur',  edge= ' lower2 ' ,  corein=9  $ 
'edgecur',  edge= ' upper2 ' ,  corein=22  $ 
'transur',  coreout=41  $ 

'edgecur',  edge= ' lowerl ' ,  corein=20  $ 
'edgecur',  edge= ' upperl ' ,  corein=19  $ 
'edgecur',  edge= ' lower2 ' ,  corein=12  $ 
'edgecur',  edge= ' upper2 ' ,  corein=21  $ 
'transur',  coreout=42  $ 

'edgecur',  edge= ' lowerl ' ,  corein=19  $ 
'edgecur',  edge= ' upperl ' ,  corein=18  $ 
'edgecur',  edge= ' lower2 ' ,  corein=ll  $ 
'edgecur',  edge= ' upper2 ' ,  corein=24  $ 
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$  'transur',  coreout=43  $ 

$  'edgecur',  edge= ' lowerl ' ,  corein=18  $ 

$  'edgecur',  edge= 'upperl ' ,  corein=17  $ 

$  'edgecur',  edge= ' lower2 ' ,  corein=10  $ 

$  'edgecur',  edge= ' upper2 ' ,  corein=23  $ 

$  'transur',  coreout=44  $ 

$  'current',  corein=41  $ 

$  'insert',  corein=42  $ 

$  'insert',  corein=43  $ 

$  'insert',  corein=44,  coreout=54  $ 

$  'getend',  corein=5,  points ' last ' ,  end=' first'  $ 

$  'line',  points=-l,  r2=l,  space=-l,  coreout=l  $ 

$  'getend',  corein=7,  point= ' last ' ,  end=' first'  $ 

$  'line',  points=-l,  r2=3,  space=-l,  coreout=3  $ 

$  'getend',  corein=8,  point= ' last ' ,  end=' first'  $ 

$  'line',  points=-l,  r2=4,  space=-l,  coreout=4  $ 

$  'edgecur',  edge= ' lowerl ' ,  corein=2  $ 

$  'edgecur',  edge= ' upperl ' ,  corein=l  $ 

$  'edgecur',  edge= ' lower2 ' ,  corein=5  $ 

$  'edgecur',  edge= 'upper2 '  ,  corein==9  $ 

$  'transur',  coreout=41  $ 

$  'edgecur',  edge= ' lowerl ' ,  corein=l  $ 

$  'edgecur',  edge= ' upperl ' ,  corein=4  $ 

$  'edgecur',  edge= ' lower2 ' ,  corein=8  $ 

$  'edgecur',  edge= ' upper2 ' ,  corein=12  $ 

$  'transur',  coreout=42  $ 

$  'edgecur',  edge= ' lowerl ' ,  corein=4  $ 

$  'edgecur',  edge= ' upperl ' ,  corein=3  $ 

$  'edgecur',  edge= ' lower2 ' ,  corein=7  $ 

$  'edgecur',  edge= ' upper2 ' ,  corein=li  $ 

$  'transur',  coreout=43  $ 

$  'edgecur',  edge= ' lowerl ' ,  corein=3  $ 

$  'edgecur',  edge= ' upperl ' ,  corein=2  $ 

$  'edgecur',  edge= ' lower2 ' ,  corein=6  $ 

$  'edgecur',  edge= ' upper2 ' ,  corein=10  $ 

$  'transur',  coreout=44  $ 

$  'current',  corein=41  $ 

$  'insert',  corein=42  $ 

$  'insert',  corein=43  $ 

$  'insert',  corein=44,  coreout=55  $ 

$  'trans',  corein=55,  origin=0. , 0. , 120. ,  coreout=56  $ 

$  'trans',  corein=51 , -56 ,  origin=250 . , 200 . , 0 . , 
coreout=51 , -56  $ 

$  'combine',  content= ' yes ' ,  corein=51 , -56 ,  fileout=l  $ 
$  'combine',  head='yes',  form='e',  triad='yes', 
corein=51 , -56 ,  fileout=2  $ 

$  'combine',  f orm= ' plot3d ' ,  corein=51 , -56 ,  fileout=3, 
f ilnam= ' aq3 . so '  $ 

$  ' end '  $ 


Figure  79  shows  the  line  segment  numbers. 
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The  very  versatile  thing  about  the  approach  taken  is  that 
once  the  first  piece  is  configured,  the  second  piece  of  its 
type  can  be  created  by  changing  the  coordinates  and  a  few 
parameters  set  at  the  beginning  of  the  run-stream. 

Bottom  surface  55  could  have  been  more  efficiently 
generated  by  first  combining  the  circular  segments  into  one 
line  segment  (lower  #2) ,  combining  the  straight  line 
segments  into  one  segment  (upper  #2) ,  and  taking  one  segment 
(say  line  segment  #1)  as  lower  #1  and  upper  #1  in  a  single 
transfinite  surface  operation.  The  analogous  thing  is,  in 
fact,  what  is  done  for  the  grid  generation  code  as  surface 
51  is  the  same  as  surface  52,  which  is  the  cut. 

The  surface  generation  data  for  the  second  well  and  the 
two  plugs  are  now  given. 


Well  0  Grid  No.  2 

$  'point',  point=l,  r=400 . , -250 . , 0 .  $ 

$  'point',  point=2,  r=400 . , 250 . , 0 .  $ 

$  'point',  point=3 ,  r=-175 . , 250. , 0.  $ 

$  'point',  point=4 ,  r=-175 . , -250 . , 0 .  $ 

$  'point',  point=5,  r=400 . , -250 . , 120 .  $ 

$  'point',  point=6,  r=400 . , 250 . , 120  .  $ 

$  'point',  point=7 ,  r=-17 5 , , 250 .  ,  120  .  $ 

$  'point',  point=8,  r=-175 . , -250 . , 120 .  $ 

$  'setnum',  segment=l,  points=13  $ 

$  'setnum',  segment=2,  points-6  $ 

$  'setnum',  segment=3,  points=25  $ 

$  'setnum',  segment=4 ,  points=17  $ 

$  'setnum',  segment=5,  points=9  $ 

$  'setval',  number=l,  value=.002  $ 

$  'setval',  number=2 ,  value-. 02  $ 

$  'setval',  number=3,  value=.04  $ 

$  'conicur',  points=-2,  type= ' circle ' ,  radius=2 . , 
angle=32 . , -32 ,  coreout=5  $ 

$  'conicur',  points=-2 ,  type= ' circle ' ,  radius=2., 
angle=125 . , 32 . ,  coreout=6  $ 


$  'conicur',  points=-2 ,  type= ' circle ' ,  radius=2 . , 
angle=235 . , 125 . ,  coreout=7  $ 

$  'conicur',  points=-2 ,  type= ' circle ' ,  radius=2 . , 
angle=-32 . , -125 . ,  coreout=8  $ 

$  'trans',  corein=5,  origin=0. , 0. , 40 , ,  coreout=49  $ 

$  'trans',  corein=5,  origin=0 . , 0 . , 120 . ,  coreout=45  $ 

$  'getend',  corein=49,  point= ' first ' ,  end=' first'  $ 

$  'getend',  corein=45,  point= ' first ' ,  end='last'  $ 

$  'line',  points=-4,  space=-2,  coreout=13  $ 

$  'getend',  corein=49,  point= ' first ' ,  end='first'  $ 

$  'getend',  corein=5,  point= ' first ' ,  end='last'  $ 

$  'line',  points=-5,  space=-3  $ 

$  'switch',  reorder= ' reversel '  $ 

$  'insert',  corein=13,  coreout=13  $ 

$  'line',  points=-3,  rl=2,  r2=6,  coreout=17  $ 

$  'getend',  corein=5,  points ' first ' ,  end=' first'  $ 

$  'line',  points=-l,  r2=2,  space=-l,  coreout=2  $ 

$  'trans',  corein=2,  origin=0 . , 0 . , 120 . ,  coreout-26  $ 

$  'edgecur',  edge= ' lowerl ' ,  corein=13  $ 

$  'edgecur',  edge- ' upperl ' ,  corein=17  $ 

$  'edgecur',  edge= ' lower2 ' ,  corein=2  $ 

$  'edgecur',  edge= ' upper2 ' ,  corein=26  $ 

$  'transur',  coreout=51  $ 

$  'current',  corein=51,  coreout=52  $ 

$  'bouncur',  corein=13  $ 

$  'rotate',  angpts=-2,  angle=0 . , -64 . ,  axcos=0 . , 0 . , 1 .  $ 

$  'switch',  reorder= ' switch ' ,  coreout=41  $ 

$  'bouncur',  corein=13  $ 

$  'rotate',  angpts=-2 ,  angle=-64 . , -157 . ,  axcos=0 . , 0 . , 1 .  $ 

$  'switch',  reorder= ' switch ' ,  coreout=42  $ 

$  'bouncur',  corein=13  $ 

$  'rotate',  angpts=-2 ,  angle=-157 . , -267 . ,  axcos=0 . , 0 . , 1 .  $ 
$  'switch',  reorder^ ' switch ' ,  coreout=43  $ 

$  'bouncur',  corein=13  $ 

$  'rotate',  angpts=-2,  angle=93 . , 0 . ,  axcos=0 . , 0 . , 1 ,  $ 

$  'switch',  reorder= ' switch ' ,  coreout=44  $ 

$  'current',  corein=41  $ 

$  'insert',  corein=42  $ 

$  'insert',  corein=43  $ 


$ 

'insert',  corein=44 

,  coreout=52 

1  $ 

$ 

'line',  points=-3. 

rl  =  l. 

r2  =  5. 

coreout=2  0 

$ 

$ 

'line',  points=-3, 

rl  =  2. 

r2  =  6. 

coreout=17 

$ 

$ 

'line',  points=-3, 

rl  =  3. 

r2=7 , 

coreout=18 

$ 

$ 

'line',  points=-3. 

11 

r2=8. 

coreout=19 

$ 

$ 

'line',  points=-2. 

rl=2 , 

r2  =  l , 

coreout=9 

$ 

$ 

'line',  points=-2, 

rl  =  3. 

r2=2 , 

coreout=10 

$ 

$ 

'line',  points=-2. 

rl=4  , 

r2  =  3. 

coreout=ll 

$ 

$ 

'line',  points=-2, 

rl  =  l. 

r2  =  4  , 

coreout=12 

$ 

$ 

'line',  points=-2. 

r  1-6 , 

r2  =  5. 

coreout=22 

$ 

$ 

'line',  points=-2 , 

rl=7. 

r2=6. 

coreout=2  3 

$ 

$ 

'line',  points=-2. 

rl=8. 

r2=7. 

coreout=24 

$ 

$ 

'line',  points=-2. 

rl  =  5. 

r2=8, 
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coreout=2 1 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 


edgecur',  edge= ' lowerl ' ,  corein 
edgecur',  edge= ' upperl ' ,  corein 
edgecur',  edge= ' lower2 ' ,  corein 
edgecur',  edge= ' upper2 ' ,  corein 
transur',  coreout=41  $ 
edgecur',  edge= ' lowerl ' ,  corein 
edgecur',  edge= ' upperl ' ,  corein 
edgecur',  edge= ' lower2 ' ,  corein 
edgecur',  edge= ' upper2 ' ,  corein 
transur',  coreout=42  $ 
edgecur',  edge= ' lowerl ' ,  corein 
edgecur',  edge= ' upperl ' ,  corein 
edgecur',  edge= ' lower2 ' ,  corein 
edgecur',  edge= 'upper2 ' ,  corein 
transur',  coreout=43  $ 
edgecur',  edge= ' lowerl ' ,  corein 
edgecur',  edge= ' upperl ' ,  corein 
edgecur',  edge= ' lower2 ' ,  corein 
edgecur',  edge==  '  upper2  '  ,  corein 
transur',  coreout=44  $ 
current',  corein=41  $ 
insert',  corein=42  $ 
insert',  corein=43  $ 
insert',  corein=44,  coreout=54 
getend ' ,  corein=5,  point='last' 
line',  points=-l,  r2=l,  space=- 
getend',  corein=7,  point='last' 
line',  points=-l,  r2=3,  space=- 
getend',  corein=8,  point='last' 
line',  points=-l,  r2=4,  space=- 


corein= 

corein= 

corein= 

corein= 


corein- 

corein= 

corein= 

corein= 


corein= 

corein= 

corein= 

corein= 


corein= 

corein= 

corein= 

corein= 


17  $ 
20  $ 

9  $ 
22  $ 

20  $ 
19  $ 
12  $ 
21  $ 

19  $ 

18  $ 
11  $ 
24  $ 

18  $ 
17  $ 

10  $ 
23  $ 


edgecur ' 
edgecur ' 
edgecur ' 
edgecur ' 
transur ' 
edgecur ' 
edgecur ' 
edgecur ' 
edgecur ' 
transur ' 
edgecur ' 
edgecur ' 
edgecur ' 
edgecur ' 
transur ' 
edgecur ' 
edgecur ' 
edgecur ' 
edgecur ' 
transur ' 
current ' 
insert ' , 


edge= ' lowerl ' 
edge= ' upperl ' 
edge= ' lower2 ' 
edge= ' upper2 ' 
coreout=41  $ 
edge= ' lowerl ' 
edge= ' upperl ' 
edge= ' lower2 ' 
edge= ' upper2 ' 
coreout=42  $ 
edge= ' lowerl ' 
edge= ' upperl ' 
edge= ' lower2 ' 
edge= ' upper2 ' 
coreout=43  $ 
edge= ' lowerl ' 
edge= ' upperl ' 
edge= ' lower2 ' 
edge= ' upper2 ' 
coreout=44  $ 
corein=41  $ 
corein=42  $ 


corein 

corein 

corein 

corein 


corein 

corein 

corein 

corein 


,  end=' first'  $ 
1,  coreout=l  $ 

,  end=' first'  $ 
1,  coreout=3  $ 

,  end=' first'  $ 
1,  coreout=4  $ 
=2  $ 

=1  $ 

=5  $ 

=9  $ 

=1  $ 

=4  $ 

=8  $ 

=  12  $ 


corein 

corein 

corein 

corein 

corein 

corein 

corein 

corein 


=  4  $ 

=  3  $ 
=7  $ 
-11  $ 

=  3  $ 
-2  $ 

=  6  $ 

=  10  $ 
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$  'insert',  corein=43  $ 

$  'insert',  corein=44,  coreout=55  $ 

$  'trans',  corein=55,  origin=0. , 0. , 120. ,  coreout=56  $ 

$  'trans',  corein=51 , -56 ,  origin=600 . , 250 . , 0 . , 
coreout=51 , -56  $ 

$  'combine',  content^ ' yes ' ,  corein=5l , -56 ,  fileout^l  $ 
$  'combine',  head='yes',  form='e',  triad='yes', 
corein=51 , -56 ,  fileout=2  $ 

$  'combine',  form== 'plot3d  '  ,  corein=51 , -56  ,  fileout=3, 
f ilnam= ' aq4 . so '  $ 

$  'end'  $ 


Plug  No.  1 


$ 

' point ' , 

point=l , 

r=l . , 0 . , 0 .  $ 

$ 

' point ' , 

point=2 , 

r=0 . , 1 . , 0 .  $ 

$ 

'point ' , 

point=3 , 

r=-l . , 0 . , 0 .  $ 

$ 

' point ' , 

point=4 , 

r=0 . , -1 . , 0 .  $ 

$ 

' point ' , 

point=5 , 

r=l . , 0 . , 80 .  $ 

$ 

' point ' , 

point=6 , 

r=0.,l.,80.  $ 

$ 

' point ' , 

point=7 , 

II 

1 

i-> 

o 

00 

o 

« 

$ 

' point ' , 

point=8 , 

• 

o 

00 

<-l 

1 

• 

o 

II 

$  'setnum',  segment=l,  points=6  $ 

$  'setnum',  segment=3 ,  points=17  $ 

$  'setval',  number=l,  value=.02  $ 

$  'conicur',  points--!,  type= ' circle ' ,  radius=l., 
angle=-141 . 3 , -48 . 8,  coreout=l  $ 

$  'conicur',  points=-l,  type= ' circle ' ,  radius=l., 
angle=-4r , 8 , 59 . 7 ,  coreout=2  $ 

$  'conicur',  points=-l,  type= ' circle ' ,  radius=l., 
angle=129 . 8 , 59 . 7 ,  coreout=3  $ 

$  'conicur',  points=-l,  type= ' circle ' ,  radius=l., 
angle=218 . 7 , 129 . 8 ,  coreout=4  $ 

$  'line',  points=-3,  rl=5,  r2=l,  space=-l  $ 

$  'switch',  reorder= ' reversel ' ,  coreout=5  $ 

$  'bouncur',  corein=5  $ 

$  'rotate',  angpts=-l,  angle=-14 1 . 3 , -48 . 8 ,  axcos=0 . , 0 . , 1 .  $ 
$  'switch',  reorder= ' switch ' ,  coreout=23  $ 

$  'bouncur',  corein=5  $ 

$  'rotate',  angpts=-l,  angle=-48 . 8 , 59 . 7 ,  axcos^O . , 0 . , 1 .  $ 

$  'switch',  reorder= ' switch ' ,  coreout=22  $ 

$  'bouncur',  corein=5  $ 

$  'rotate',  angpts=-l,  angle=129 . 8 , 59 . 7 ,  axcos=0 . , 0 . , 1 .  $ 

$  'switch',  reorder^ ' switch ' ,  coreout=24  $ 

$  'bouncur',  corein=5  $ 

$  'rotate',  angpts=-l,  angle=2 18 . 7 , 129 . 8 ,  axcos=0 . , 0 . , 1 .  $ 

$  'switch',  reorder= ' switch ' ,  coreout=21  $ 

$  'edgecur',  edge= ' lowerl ' ,  corein=4  $ 

$  'edgecur',  edge= ' upperl ' ,  corein=2  $ 

$  'edgecur',  edge= ' lower2 ' ,  corein^l  $ 

$  'edgecur',  edge= ' upper2 ' ,  corein=3  $ 
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$  'transur',  coreout=25  $ 

$  'trans',  origin=0 . , 0 . , 80 . ,  coreout=26  $ 

$  'trans',  corein=21 , -26 ,  origin=250. , 200. , 0. , 
coreout=21, -26  $ 

$  'combine',  content= ' yes ' ,  corein=21 , -26 ,  fileout=l  $ 
$  'combine',  head='yes',  form='e',  triad='yes', 
corein=21 , -26 ,  fileout=2  $ 

$  'combine',  form='plot3d' ,  corein=21, -26,  fileout=3, 
filnam^'plugl.so'  $ 

$  'end'  $ 


Plug  No.  2 


$ 

' point ' , 

point^l , 

r=2 . , 0. 

$ 

' point ' , 

point=2 , 

r=0 . , 2 . 

$ 

' point ' , 

point=3 , 

o 

fM 

1 

II 

u 

$ 

' point ' , 

point=4 , 

r-0 . , -2 

$ 

' point ' , 

point=5 , 

r=2 . , 0 . 

$ 

' point ' , 

point=6 , 

r=0 . , 2 . 

$ 

'point ' , 

point=7 , 

r=-2 . , 0 

$ 

' point ' , 

point=8 , 

r=0. , -2 

' setnum ' , 
' setnum ' , 
' setval ' , 
' conicur ' 


conicur ' 


)oint=l,  r=2.,0.,0.  $ 

)oint=2,  r=0.,2.,0.  $ 

)oint=3,  r=-2.,0.,0,  $ 

)oint=4 ,  r=0.,-2.,0.  $ 

)oint=5,  r=2.,0.,40.  $ 

)oint=6,  r=0.,2.,40.  $ 

5oint=7 ,  r=-2.,0.,40.  $ 

3oint=8,  r=0.,-2.,40.  $ 
segment=l,  points-6  $ 
segment=3,  points=9  $ 
number=l,  value=.04  $ 
points=-l,  type= ' circle ' ,  radius=2., 
angle=-125 . , -32 . ,  coreout=l  $ 
points=-l,  type= ' circle ' ,  radius=2.. 


' conicur ' 


,32.,  coreout=2 
type= ' circle ' , 
,32.,  coreout=3 
type= 'circle ' , 


radius=2 

$ 

radius=2 


,125.,  coreout=4  $ 


angle=-32 . , 32 . ,  coreout=2  $ 

'conicur',  points=-l,  type= ' circle ' ,  radius=2., 
angle=125 . , 32 . ,  coreout=3  $ 

'conicur',  points=-l,  type= ' circle ' ,  radius=2 . , 
angle=235 . , 125 . ,  coreout=4  $ 

'line',  points=-3,  rl=5,  r2=l,  space=-l  $ 

'switch',  reorder= ' reversel ' ,  coreout=5  $ 

'bouncur',  corein=5  $ 

'rotate',  angpts=-l,  angle=-125 . , -32 . ,  axcos=0 . , 0 . , 1 . 
'switch',  reorder^ ' switch ' ,  coreout=23  $ 

'bouncur',  corein=5  $ 

'rotate',  angpts=-l,  angle=-32 . , 32 . ,  axcos=0 . , 0 . , 1 .  $ 
'switch',  reorder= ' switch ' ,  coreout=22  $ 

'bouncur',  corein=5  $ 

'rotate',  angpts=-l,  angle=125 . , 32 . ,  axcos=0 . , 0 . , 1 .  $ 
'switch',  reorder= ' switch ' ,  coreout=24  $ 

'bouncur',  corein=5  $ 

'rotate',  angpts=-l,  angle=235 . , 125 . ,  axcos=0 . , 0 . , 1 . 
'switch',  reorder= ' switch ' ,  coreout=21  $ 


' bouncur ' 
' rotate ' , 
' switch ' , 
' bouncur ' 
' rotate ' , 
' switch ' , 
' bouncur ' 
' rotate ' , 
' switch ' , 
' bouncur ' 
' rotate ' , 
' switch ' , 
' edgecur ' 
' edgecur ' 
' edgecur ' 
' edgecur ' 
' transur ' 


edge= ' lowerl ' , 
edge= ' upperl ' , 
edge= ' lower2 ' , 
edge= ' upper2 ' , 
coreout=25  $ 


corein=4 

corein=2 

corein=l 

corein=3 
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$  'trans',  origin=0 . , 0 . , 40 . ,  coreout=26  $ 

$  'trans',  corein=21 , -26 ,  origin=600. , 250. , 0. , 
coreout=21 , -26  $ 

$  'combine',  content= ' yes ' ,  corein=2 1 , -26 ,  fileout=l  $ 
$  'combine',  head='yes',  form='e’,  triad='yes', 
corein=21, -26,  fileout=2  $ 

$  'combine',  form= 'plot3d ' ,  corein=21 , -26 ,  fileout=3, 
f ilnam= ' plug2 . so '  $ 

$  'end'  $ 


Grid  Generation  Data 


The  data  for  the  portion  of  EAGLE  that  generates  the 
grid  from  the  produced  surfaces  is  the  same  for  each  piece, 
except  for  the  numbers  on  the  faces  and  the  output  file 
names.  This  data  for  the  first  O  grid  is  now  given. 


$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 


'  initial 

'  file'  , 

'  block ' $ 
'point' , 
'  point ' , 
'point ' , 
'  point ' , 

'  point ' , 
'  point ' , 
'  point ' , 
'point ' , 


',  kstore= ' core ' ,  accpar= ' optimum ' ,  itmax=100, 
tol=. 000001,  contyp= ' initial '  $ 


file=ll  , 

all= 'yes ' 

$ 

point=l. 

locat=l , 1 

,1  $ 

point=2 , 

opoint=l , 

segment=55 , 

direct=l , 

ndex=l 

point=3 , 

opoint=2 , 

segment=55 , 

direct=2 , 

ndex=2 

point=4 , 

opoint=3 , 

segment=55 , 

direct=-l 

f 

ndex=l  $ 
point=5 , 

opoint=l , 

segment=53 , 

direct=3 , 

ndex=2 

point=6 , 

opoint=5 , 

segment=56 , 

direct=l , 

ndex=l 

point=7 , 

opoint=6. 

segment=56 , 

direct=2 , 

ndex=2 

point=8 , 
ndex=l  $ 

opoint=7 , 

segment=56 , 

direct=-l 

f 

$  ' size ' ,  size=7  $ 


$ 

'segment',  segment=51. 

start=l , 

end=8 , 

class= 'fix' 

$ 

$ 

' segment ' ,  segment=52 , 

start=2 , 

end=7 , 

class= 'fix' 

$ 

$ 

'segment',  segment=53. 

start=l , 

end=6 , 

class= 'fix' 

$ 

$ 

'segment',  segment=54. 

start=4 , 

end=7 , 

class= 'fix' 

$ 

$ 

' segment ' ,  segment=55 , 

start=l , 

end=3 , 

class= 'fix' 

$ 

$ 

'segment',  segment=56. 

start=5 , 

end=7 , 

class= 'fix' 

$ 

$  'store',  file=14,  outer= ' plot3d ' ,  f ilnam= ' aq3 . egl '  $ 
$  'end'  $ 

$  ' error '$ 

$  'end'  $ 


$ 

$ 


$ 

$ 

$ 


The  point  and  segment  (face)  numbers  are  shown  in  the  topo¬ 
logical  configuration  given  in  Figure  80. 
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Figure  80.  Topology  for  Grid  Program 

FEM  Grid  Preparation  Data 

The  question-and-answer  sequence  for  the  program  to 
combine  the  four  separate  grids  into  the  proper  FEM  format 
is  now  given. 

COMMAND? 

combine 

NUMBER  OF  FILES 
4 

FILE  NAME  NO.  1  ? 
aq3 . egl 

FOR  BLOCK  1  MATERIAL  TYPE  NUMBER,  THETA,  PHI,  AND  PSI? 

10  0  0 
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REMOVING 

DUPLICATE 

NODES, 

N 

= 

1000 

NUMNP  = 

0 

REMOVING 

DUPLICATE 

NODES , 

N 

= 

2000 

NUMNP  = 

0 

REMOVING 

DUPLICATE 

NODES , 

N 

= 

3000 

NUMNP  = 

0 

REMOVING 

DUPLICATE 

NODES , 

N 

= 

4000 

NUMNP  = 

0 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

5000 

NUMNP  = 

0 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

6000 

NUMNP  = 

0 

FILE  NAME  NO.  2  ? 
q4 . egl 

FOR  BLOCK  1  MATERIAL  TYPE  NUMBER,  THETA,  PHI,  AND  PSI? 
10  0  0 


REMOVING 

DUPLICATE 

NODES, 

N 

= 

1000 

NUMNP  = 

6500 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

2000 

NUMNP  = 

6500 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

3000 

NUMNP  = 

6500 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

4000 

NUMNP  = 

6500 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

5000 

NUMNP  = 

6500 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

6000 

NUMNP  = 

6500 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

7000 

NUMNP  = 

6500 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

8000 

NUMNP  = 

6500 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

9000 

NUMNP  = 

6500 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

10000 

NUMNP  = 

6500 

REMOVING 

DUPLICATE 

NODES , 

N 

= 

11000 

NUMNP  = 

6500 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

12000 

NUMNP  = 

6500 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

13000 

NUMNP  = 

6500 

FILE  NAME  NO.  3  ? 
plugs . egl 

FOR  BLOCK  1  MATERIAL  TYPE  NUMBER,  THETA,  PHI,  AND  PSI? 
10  0  0 


REMOVING 

DUPLICATE 

NODES , 

N 

1000 

NUMNP 

= 

12850 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

2000 

NUMNP 

= 

12850 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

3000 

NUMNP 

= 

12850 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

4000 

NUMNP 

= 

12850 

REMOVING 

DUPLICATE 

NODES , 

N 

= 

5000 

NUMNP 

= 

12850 

REMOVING 

DUPLICATE 

NODES , 

N 

6000 

NUMNP 

= 

12850 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

7000 

NUMNP 

= 

12850 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

8000 

NUMNP 

= 

12850 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

9000 

NUMNP 

= 

12850 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

10000 

NUMNP 

= 

12850 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

11000 

NUMNP 

= 

12850 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

12000 

NUMNP 

= 

12850 

REMOVING 

DUPLICATE 

NODES, 

N 

13000 

NUMNP 

= 

12850 

FILE  NAME  NO.  4  ? 
plug4 . egl 

FOR  BLOCK  1  MATERIAL  TYPE  NUMBER,  THETA,  PHI,  AND  PSI? 
10  0  0 
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REMOVING 

DUPLICATE 

NODES, 

N 

= 

1000 

NUMNP  = 

13122 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

2000 

NUMNP  = 

13122 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

3000 

NUMNP  = 

13122 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

4000 

NUMNP  = 

13122 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

5000 

NUMNP  = 

13122 

REMOVING 

DUPLICATE 

NODES, 

N 

=r 

6000 

NUMNP  = 

13122 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

7000 

NUMNP  = 

13122 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

8000 

NUMNP  = 

13122 

REMOVING 

DUPLICATE 

NODES, 

N 

9000 

NUMNP  = 

13122 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

10000 

NUMNP  = 

13122 

REMOVING 

DUPLICATE 

NODES , 

N 

= 

11000 

NUMNP  = 

13122 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

12000 

NUMNP  = 

13122 

REMOVING 

DUPLICATE 

NODES, 

N 

= 

13000 

NUMNP  = 

13122 

FINAL  NUMBER  ON  NODES  =  13266 


COMMAND? 

bandwidth 


INITIAL  MAXIMUM  DIFFERENCE  =  12881 


ITERATION 

= 

1 

BEST 

MAXIMUM 

DIFFERENCE 

= 

1436 

ITERATION 

7 

BEST 

MAXIMUM 

DIFFERENCE 

1413 

ITERATION 

= 

8 

BEST 

MAXIMUM 

DIFFERENCE 

= 

1375 

ITERATION 

= 

9 

BEST 

MAXIMUM 

DIFFERENCE 

= 

1331 

ITERATION 

= 

10 

BEST 

MAXIMUM 

DIFFERENCE 

= 

1288 

ITERATION 

= 

11 

BEST 

MA.XIMUM 

DIFFERENCE 

= 

1252 

ITERATION 

= 

12 

BEST 

MAXIMUM 

DIFFERENCE 

= 

1222 

ITERATION 

13 

BEST 

MAXIMUM 

DIFFERENCE 

= 

1220 

COMMAND? 

be 

BLOCK  NUMBER,  IIB,  JIB,  KIB,  I2B,  J2B,  K2B,  IBC,  BV? 

1  6  13  1  11  13  25  1  0 

COMMAND? 

be 

BLOCK  NUMBER,  IIB,  JIB,  KIB,  I2B,  J2B,  K2B,  IBC,  BV? 

2  6  13  1  11  13  25  1  0 

COMMAND? 

be 

BLOCK  NUMBER,  IIB,  JIB,  KIB,  I2B,  J2B,  K2B,  IBC,  BV? 
1  1  1  17  21  1  25  -2  -.4 

COMMAND? 

be 
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BLOCK  NUMBER,  IIB,  JIB,  KIB,  I2B,  J2B,  K2B,  IBC,  BV? 
2  1  1  9  21  1  25  -2  -.2 


COMMAND? 

out 


OUTPUT  FILE  NAME  FOR  SEEPAGE/GROUNDWATER  MODEL  (WITHOUT 

EXTENSIONS) ? 

aqbs 

TITLE  CARD? 

Two  Wells  in  an  Aquifer 

DATUM? 

120 

Kl,  K2,  AND  K3  FOR  MATERIAL  TYPE  1  ? 

.1  .1  .1 

COMMAND? 

end 


FEM  Grid 

Portions  of  the  FEM  grid  data  are  now  given. 


Two  Wells  in  an  Aquifer 


13266 

12120 

1  480 

120.00 

1 

0. 

10  0. 

10  0. 

10 

1 

0 

249.00 

199.95 

0.00 

0.00 

2 

0 

248.50 

199.95 

0.00 

0.00 

3 

0 

249.07 

199.64 

0.00 

0.00 

4 

0 

249.00 

199.95 

8.07 

0.00 

5 

0 

249.03 

200.25 

0.00 

0 . 00 

6 

0 

249 . 36 

199.87 

0.00 

0.00 

7 

0 

248.57 

199.44 

0.00 

0.00 

8 

0 

248.36 

199.94 

8.11 

0.00 

9 

0 

248.53 

200.45 

0.00 

0.00 

10 

0 

247.25 

199.95 

0.00 

0.00 

13261 

0 

602.70 

250.97 

111.93 

0.00 

13262 

-1 

601.89 

250.66 

120.00 

0.00 

13263 

0 

602.78 

250.32 

120.00 

0.00 

13264 

0 

604 . 78 

250.57 

120.00 

0.00 

13265 

-1 

601.70 

251.06 

120.00 

0.00 

13266 

0 

602 . 68 

250.96 

120.00 

0. 00 
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1 


1 

394  287 

395 

469 

493 

378 

505 

568 

0.00 

0.00 

0.00 

2 

287  200 

288 

395 

378 

278 

385 

505 

0.00 

0.00 

0.00 

3 

200  132 

201 

288 

278 

197 

285 

385 

0.00 

0.00 

0.00 

4 

132  179 

263 

201 

197 

262 

366 

285 

0.00 

0.00 

0.00 

5 

179  251 

290 

263 

262 

354 

398 

366 

0.00 

0.00 

0.00 

6 

251  169 

202 

290 

354 

252 

291 

398 

0.00 

0.00 

0.00 

7 

169  106 

133 

202 

252 

170 

203 

291 

0.00 

0.00 

0.00 

8 

106  60 

81 

133 

170 

107 

134 

203 

0.00 

0.00 

0.00 

9 

60  29 

44 

81 

107 

61 

82 

134 

0.00 

0.00 

0.00 

10 

29  11 

20 

44 

61 

30 

45 

82 

0.00 

0.00 

. 

0.00 

12113 

9794  10100 

10371  10080 

10099 

10388 

10646 

10370 

0.00 

0.00 

0.00 

12114 

10100  10389 

10647  10371 

10388 

10662 

10909 

10646 

0.00 

0.00 

0.00 

12115 

10389  10663 

10910  10647 

10662 

10923 

11163 

10909 

0.00 

0.00 

0.00 

12116 

9447  9772 

10010  ' 

9687 

9771 

10079 

10309 

10003 

0.00 

0.00 

0.00 

12117 

9772  10080 

10315  10010 

10079 

10370 

10598 

10309 

0.00 

0.00 

0.00 

12118 

10080  10371 

10603  10315 

10370 

10646 

10872 

10598 

0 . 00 

0.00 

0.00 

12119 

10371  10647 

10876  10603 

10646 

10909 

11134 

10872 

0.00 

0.00 

0.00 

12120 

10647  10910 

11137  10876 

10909 

11163 

11385 

11134 

0.  0 

• 

0. 

4168 

3860  4176 

4500 

-0 

.40 

3860 

3576  3878 

4176 

-0 

.40 

3576 

3303  3594 

3878 

-0 

.40 

3303 

3535  3836 

3594 

-0 

.40 

3535 

3828  4142 

3836 

-0 

.40 

3828 

3529  3829 

4142 

-0 

.40 

3529 

3240  3530 

3829 

-0 

.40 

3240 

2958  3241 

3530 

-0 

.40 

12885 

12775  12895 

12988 

-0 

.20 

12775 

12892  12993 

12895 

-0 

.20 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

1 

1 

1 

1 
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12892  12991  13073 
12991  13045  13110 
13045  13117  13166 
13117  13172  13208 
13172  13213  13238 
13213  13242  13257 
13242  13260  13265 


12993 

-0.20 

13073 

-0.20 

13110 

-0.20 

13166 

-0.20 

13208 

-0.20 

13238 

-0.20 

13257 

-0.20 
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WATERWAYS  EXPERIMENT  STATION  REPORTS 
PUBLISHED  UNDER  THE  COMPUTER-AIDED 
STRUCTURAL  ENGINEERING  (CASE)  PROJECT 


I  (  Clinical  Report  K-70-1 
Irr  Unction  Report  0-70 

Technical  Report  K-80-1 
Ttichnical  Report  ‘^-00-2 

Instruction  Report  K-00-1 

Instruction  Report  K-80-3 
Instruction  Report  K-80-4 

Instruction  Report  K-80'6 
Instruction  Report  K-80-7 
Technical  Report  K-80-4 

Technical  Report  K-80-5 
Instruction  Report  K-81-2 

Instruction  Report  K-81-3 
Instruction  Report  K-81-4 
Instruction  Report  K-81-6 

Instruction  Report  K-81-7 
Instruction  Report  K-81-9 
Technical  Report  K-81-2 
Instruction  Report  K-82-6 


Title 


Oatft 


List  of  Computer  Programs  for  Computer-Aided  Structural  Engineering 

User's  Guido;  Computer  Program  with  Interactive  Graphics  for 
Analysis  of  Plane  Frame  Structures  (CFRAME) 

Suivoy  of  Bridge-Oriented  Design  Software 

Evaluation  of  Computer  Programs  for  the  Design/Analysis  of 
Highway  and  Railway  Bridges 

User's  Guido:  Computer  Program  for  Design/Review  of  Curvi¬ 
linear  Conduits/Culverts  (CURCON) 

A  Three-Dimensional  Finite  Element  Data  Edit  Program 

A  Throe-Dimensional  Stability  Analysis/Design  Program  (3DSAD) 
Report  1 :  General  Geometry  Module 
Report  3:  General  Analysis  Module  (CGAM) 

Report  4:  Special-Purpose  Modules  for  Dams  (CDAMS) 

Basic  User's  Guide:  Computer  Program  for'Design  and  Analysis 
of  Inverted-T  Retaining  Walls  and  Floodwalls  (TWDA) 

User's  Reference  Manual:  Computer  Program  for  Design  and 
Analysis  of  Inverted-T  Retaining  Walls  and  Roodwalls  (TWDA) 

Documentation  of  Finite  Element  Analyses 
Report  1 :  Longview  Outlet  Works  Conduit 
Report  2:  Anchored  Wall  Monolith,  Bay  Springs  Lock 

Basic  Pile  Group  Behavior 

User's  Guide:  Computer  Program  for  Design  and  Analysis  of  Sheet 
Pile  Walls  by  Classical  Methods  (CSHTWAL) 

Report  1 :  Computational  Processes 
Report  2:  Interactive  Graphics  Options 

Validation  Report:  Computer  Program  for  Design  and  Analysis  of 
Inverted-T  Retaining  Walls  and  Floodwalls  (TWDA) 

User's  Guide:  Computer  Program  for  Design  and  Analysis  of 
Cast-in- Place  Tunnel  Linings  (NEWTUN) 

User's  Guide:  Computer  Program  for  Optimum  Nonlinear  Dynamic 
Design  of  Reinforced  Concrete  Slabs  Under  Blast  Loading 
(CBARCS) 

User's  Guide:  Computer  Program  for  Design  or  Investigation  of 
Orthogonal  Culverts  (CORTCUL) 

User's  Guide:  Computer  Program  for  Three-Dimensional  Analysis 
of  Building  Systems  (CTABS80) 

Theoretical  Basis  for  CTABS80:  A  Computer  Program  (or 
Three-Dimensional  Analysis  of  Building  Systems 

User's  Guide:  Computer  Program  (or  Analysis  of  Beam-Column 
Stoictures  with  Nonlinear  Supports  (CBEAMC) 


reb 

(Aar 

Jan  1380 
Jan  198C 

Feb  1980 

Mar  1980 

Jun  1980 
Jun  :S32 
Aug  1983 

Dec  1980 

Dec  1980 


Dec  1980 
Dec  1930 

Dec  1930 


Feb  1981 
Msr  1981 

Feb  1931 
Mar  1981 
Mar  1981 

Mar  1981 
Aug  1961 
Sep  1981 
JUn1982 
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WATERWAYS  EXPERIMENT  STATION  REPORTS 
PUBLISHED  UNDER  THE  COMPUTER-AIDED 
STRUCTURAL  ENGINEERING  (CASE)  PROJECT 


Instruction  Report  K-02-7 
Instruction  Report  K-03-1 
Instruction  Report  K-83-2 
Instruction  Report  K-83-5 

Technical  Report  K-83-1 
Technical  Report  K-83-3 

Technical  Report  K-83-4 
Instruction  Report  K-84-2 

Instruction  Report  K-84-7 

Instruction  Report  K-84-8 

Instruction  Report  K-84-1 1 

Technical  Report  K-84-3 

Technical  Report  ATC-86'5 

Technical  Report  ITL-87-2 
Instruction  Report  ITL-87- 1 
Instruction  Report  ITL-87-2 
Technical  Report  ITL-87-6 
Instruction  Report  ITL-87-3 


(Continued) 

Title 

User's  Guide:  Computer  Program  for  Bearing  Capacity  Analysis 
of  Shallow  Foundations  (CBEAR) 

User's  Guide:  Computer  Program  with  Interactive  Graphics  for 
Analysis  of  Plane  Frame  Structures  (CFRAME) 

User's  Guide:  Computei  Program  for  Generation  of  Engineering 
Geometry  (SKETCH) 

User's  Guide:  Computer  Program  to  Calculate  Shear,  Moment, 
and  Thrust  (CSMT)  from  Stress  Results  of  a  Two-Dimensional 
Finite  Element  Analysis 

Basic  Pile  Group  Behavior 

Reference  Manual;  Computer  Graphics  Program  for  Generation  of 
Engineering  Geometry  (SKETCH) 

Case  Study  of  Six  Major  General-Purpose  Finite  Element  Programs 

User's  Guide:  Computer  Program  for  Optimum  Dynamic  Design 
of  Nonlinear  Metal  Plates  Under  Blast  Loading  (CSDOOR) 

User's  Guide:  Computer  Program  for  Determining  Induced 
Stresses  and  Consolidation  Settlements  (CSETT) 

Seepage  Analysis  of  Confined  Flow  Problems  by  the  Method  of 
Fragments  (CFRAG) 

User's  Guide  for  Computer  Program  CGFAG,  Concrete  General 
Flexure  Analysis  with  Graphics 

Computer-Aided  Drafting  and  Design  for  Corps  Structural 
Engineers 

Decision  Logic  Table  Formulation  of  ACI  318-77,  Building  Code 
Requirements  for  Reinforced  Concrete  for  Automated  Con¬ 
straint  Processing,  Volumes  I  and  II 

A  Case  Committee  Study  of  Finite  Element  Analysis  of  Concrete 
Flat  Slabs 

User's  Guide:  Computer  Program  for  Two-Dimensional  Analysis 
of  U-Frame  Structures  (CUFRAM) 

User's  Guide:  For  Concrete  Strength  Investigation  and  Design 
(CASTR)  in  Accordance  with  ACI  318-83 

Finite-Element  Method  Package  for  Solving  Steady-State  Seepage 
Problems 

User's  Guide:  A  Three  Dimensional  Stability  Analysis/Design 
Program  (3DSAD)  Module 
Report  1 :  Revision  1 :  General  Geometry 
Report  2:  General  Loads  Module 
Report  6;  Free-Body  Module 


Date 

Jun  1982 
Jan  1983 

Jun  1983 

Jul  1983 

Sep  1983 
Sep  1983 

Oct  1983 
Jan  1984 

Aug  1984 

Sep  1984 

Sep  1984 

Oct  1984 

Jun  1986 

Jan  1987 
Apr  1987 
May  1987 
May  1987 

Jun  1987 

Jun  1907 
Sep  1989 
Sep  1989 
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Instruction  Report  ITL-87-4 
Technical  Report  ITL-37-4 


instruction  Report  GL-87-1 

Instruction  Report  ITL-87-5 
Instruction  Report  ITL-87-6 

Technical  Report  ITL-87-8 

Instruction  Report  ITL-88-1 

Technical  Report  ITL-88-1 

Technical  Report  ITL-88-2 

Instruction  Report  ITL-88-2 

Instruction  Report  ITL-88-4 

Instruction  Report  GL-87-1 

Technical  Report  ITL-89-3 
Technical  Report  ITL-89-4 


(Continued) 

Title 


User’s  Guide:  2-D  Frame  Analysis  Link  Program  (LINK2D) 


Finite  Element  Studies  of  a  Horizontally  Framed  Miter  Gate 

Report  1 :  Initial  and  Refined  Rnite  Element  Models  (Phases 
A,  B,  and  C),  Volumes  I  and  II 
Report  2:  Simplified  Frame  Model  (Phase  D) 

Report  3:  Alternate  Configuration  Miter  Gate  Finite  Element 
Studies-Open  Section 

Report  4:  Alternate  Configuration  Miter  Gate  Finite  Element 
Studies-Closed  Sections 

Report  5:  Alternate  Configuration  Miter  Gate  Finite  Element 
Studies-Additional  Closed  Sections 
Report  6;  Elastic  Buckling  of  Girders  in  Horizontally  Framed 
Miter  Gates 

Report  7:  Application  and  Summary 


User’s  Guide:  UTEXAS2  Slope-Stability  Package:  Volume  I, 
User’s  Manual 


Sliding  Stability  of  Concrete  Structures  (CSLIDE) 

Criteria  Specifications  for  and  Validation  of  a  Computer  Program 
for  the  Design  or  Investigation  of  Horizontally  Framed  Miter 
Gates  (CMITER) 

Procedure  for  Static  Analysis  of  Gravity  Dams  Using  the  Finite 
Element  Method  -  Phase  1a 

User’s  Guido:  Computer  Prograi.i  for  Analysis  of  Planar  Grid 
Structures  (CGRID) 

Development  of  Design  Formulas  for  Ribbed  Mat  Foundations 
on  Expansive  Soils 

User’s  Guide:  Pile  Group  Graphics  Display  (CPGG)  Post¬ 
processor  to  CPGA  Program 

User’s  Guide  for  Design  and  Investigation  of  Horizontally  Framed 
Miter  Gates  (CMITER) 

User’s  Guide  for  Revised  Computer  Program  to  Calculate  Shear, 
Moment,  and  Thrust  (CSMT) 

User’s  Guide:  UTEXAS2  Slope-Stability  Package:  Volume  II, 
Theory 

User’s  Guide:  Pile  Group  Analysis  (CPGA)  Computer  Group 

CBASIN-Structural  Design  of  Saint  Anthony  Falls  Stilling  Basins 
According  to  Corps  of  Engineers  Criteria  for  Hydraulic 
Structures:  Computer  Program  X0098 


Date 
Jun  1987 
Aug  1987 


Aug  1987 

Oct  1987 
Dec  1987 

Jan  1988 

Feb  1988 

Apr  1988 

Apr  1988 

Jun  1988 

Sep  1988 

Feb  1989 

Jul1989 
Aug  1989 


(Continued) 


WATERWAYS  EXPERIMENT  STATION  REPORTS 
PUBLISHED  UNDER  THE  COMPUTER-AIDED 
STRUCTURAL  ENGINEERING  (CASE)  PROJECT 


Technical  Report  ITL-89-5 

Technical  Report  ITL-89-6 
Contract  Report  ITL-89-1 
instruction  Report  ITL-90-1 
Technical  Report  ITL-90-3 

Instruction  Report  ITL-90-6 
Instruction  Report  ITL-90-2 
Technical  Report  ITL-91-3 


(Concluded) 

Title 

CCHAN-Structural  Design  of  Rectangular  Channels  According 
to  Corps  of  Engineers  Criteria  for  Hydraulic 
Structures:  Computer  Program  X0097 

The  Response-Spectrum  Dynamic  Analysis  of  Gravity  Dams  Using 
the  Finite  Element  Method;  Phase  li 

State  of  the  Art  on  Expert  Systems  Applications  in  Design, 
Construction,  and  Maintenance  of  Structures 

User's  Guide:  Computer  Program  for  Design  and  Analysis 
of  Sheet  Pile  Walls  by  Classical  Methods  (CWALSHT) 

Investigation  and  Design  of  U-Frame  Structures  Using 
Program  CUFRBC 

Volume  A:  Program  Criteria  and  Documentation 
Volume  B:  User’s  Guide  for  Basins 
Volume  C:  User's  Guide  for  Channels 

User's  Guide:  Computer  Program  for  Two-Dimensional  Analysis 
of  U-Frame  or  W-Frame  Structures  (CWFRAM) 

User’s  Guide:  Pile  Group>-Concrete  Pile  Analysis  Program 
(CPGC)  Preprocessor  to  CPGA  Program 

Application  of  Finite  Element,  Grid  Generation,  and  Scientific 
Visualization  Techniques  to  2-D  and  3-D  Seepage  and 
Groundwater  Modeling 


Date 

Aug  1989 

Aug  1 989 
Sep  1989 
Feb  1990 
May  1990 

Sep  1990 
Jun  1990 
Sep  1990 


